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ABSTRACT 


The development of a universal solution of the main problem in artificial satellite 
theory has only recently been accomplished with the aid of high powered computers. 
The solution-to this Jong standing problem is an analytical expression that is similar in 
form to the two-body solution. An analysis is presented in which the solution is com- 
pared with the two-body solution, a proven numerical solution, and actual measured 
satellite data. The solution is shown to be significantly more accurate than the two-body 
solution. The theoretical accuracy of the solution is confirmed. The solution compares 
extremely well with a proven numerical solution for at least 41 orbits with a relative error 
on the order of J’@. The solution compares extremely well with measured satellite data 
for satellites in near Earth orbits. For a satellite in orbit at an altitude of approximately 
1000-kilometers, the solution reduces the error of the two-body solution by about 95%. 
For satellites in orbit at semisynchronous and geosynchronous altitudes, the solution 
reduces the error of the two-body solution by at least 50%. The solution is free of 
singularities and is valid for all eccentricities and inclinations. 
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I. INTRODUCTION 


With the development of any new method or theory the question always arises on 
whether the approach is valid or practical for ordinary use. This is particularly true in 
the prediction of satellite motion. Ever since Sir Issac Newton’s discovery of the law 
of universal gravitation, new methods have been developed to better predict the motion 
of the heavenly bodies, Usually the method contains one or more restrictions that limits 
the practical use of the solution. The goal, of course, is to develop a solution that is 
valid and possesses no restrictions. Recently, such a solution has been formulated. 

This analysis continues the work that was begun by Snider [Ref. 1] and Sagovac 
[Ref. 2]. From their research, a higher order universal solution of the motion of an 
artificial satellite around an oblate planet was developed. The solution is free of 
singularities and is theoretical valid for all orbital parameters. The purpose and scope 
of this work is to compare the solution with proven numerical solutions and actual 
measured satellite data in order to determine whether the theoretical work is valid and 
practical. 

The first chapter summarizes the development of the theory and presents the sol- 
ution in its entirety. Also-included is a somewhat less accurate simplified solution. An 
explanation of the solution near the critical inclinations is presented. The chapter con- 
cludes with a discussion on the conservation of specific mechanical energy. The next 
chapter describes the method of analysis and explains the type of integration routine 
exercised in the evaluation. The method of comparison is presented next in which the 
error parameters are described in detail. The results follow which include both a detailed 
discussion and a graphic representation. 

The analysis is objective in nature and designed to demonstrate both the advantages 
and disadvantages of the theory and the solution. Before the solution can be applied 
extensively, a general understanding of its strengths and weaknesses must be determined. 














Ii. BACKGROUND 


A. ORBITAL KINEMATICS 


A reference system for a planet in spherical coordinates (r,«, 8) is shown in Figure 
1. The radial distance. (r)-is measured from the center of the planet (O) to the satellite 
(S). The line (O y) is in a direction fixed with respect to an inertial coordinate system. 
For Earth, the line (O y) is in the direction of the vernal equinox. The right ascension 
angle («) is measured in the planet’s equatorial plane eastward from the line (O y). The 
declination or latitude (f) is the angle measured northward from the equator. The po- 
sition vector (r) of the satellite in the spherical coordinate system is: 


r = r(cosacosf)b,+r(sinacos 8) b, +r (sin B ) b3 (1) 


where (b, , b, ,b, ) are orthonormal base vectors fixed in the direction shown. 


polor axis 






equatoriol- 
plane 


Figure 1. Spherical coordinate system. 


A reference system for a satellite in polar coordinates (r, 0) within a rotating orbital 
plane is shown in Figure 2. The satellite’s position and velocity vectors are contained 
within the orbital plane. The argument of latitude (@) is measured in the orbitai piane 


two 

















from the ascending node to the satellite. The inclination (i) of the orbital plane is the 
angle measured between the equatorial plane and the orbital plane. The longitude of the 
ascending node (Q) is measured from the line (O y) to the ascending node. The as- 
cending node lies on the line of nodes which form the intersection of the equatorial and 
orbital planes. 


orbital plane 






equatorial plane 


Figure 2, Orbital plane. 


The basis (b,,b,,b,;) may be transformed into another orthonormal basis 
( B,, B, ,B,) by a succession of three rotations. First the basis (b, , b, , b, ) is rotated 
about the b, direction by the angle ©. The basis is then rotated about the new first 
coordinate vector by the angle i. The final rotation is about the new third coordinate 
vector by the angle 0. The position vector (r) has only one component in the rotating 
basis. 


r= 7B (2) 











The components of r in the fixed’ basis are: 


r = r(cos@cosQ-—sin@cosisinQ)b, 4 
+ r(cos@ sin Q+-sin@ cos icosQ)b, +r (sin @ sini) bs (3) 


Equating the components of equations (1) and (3), the following relations among 
the angles («, #) of the spherical. coordinate system and the astronomical angles (i, 
Q,6) can be obtained. 


snf = sin@sini (4) 
cosB = cos@ sect 2 Q) (5) 


The velocity (dr/dt) of the satellite is cbtained by differentiating equation (2) with 
respect to time (¢) which results in: 
dé 


Spar @ (14 tand cot 716 )B (6) 


deo dr 
dt dt 


Equations (2) and (6) represent the orbital kinematics of a satellite in the polar co- 
ordinate system. The position and velocity vector expressions describe the motion of a 
satellite in a particular reference system and provide the information needed to develop 
the equations of motion in that system. These expressions are referenced to the zrue, 
rather than mean, orbital plane and were originally formulated by Struble 
[Ref.3,4,5]. 


B. EQUATIONS OF MOTION 


The motion of all objects is mathematically described by the equations of motion 
that govern them. For an oblate planet, the expressions for the kinetic and potential 
energies per unit mass of an orbiting satellite in spherical coordinates are respectively: 


r= t]{(£)er(B)sreor(#)] oo 
2 

-_ Sit] 14 (1-3) (8) 
} 











In the above equations, (Jf) is the mass of the planet, (G) is the universal 
gravitational constant, (R) is the equatorial radius of the planet, and (J,) is the coef- 
< ficient of the second zonal harmonic of the planet’s gravitational potential. The gov- 
erning equations of motion can be determined by substituting equations (7) and (8) into 
Lagrange’s equations which are represented by: 


Vii We Gla 6 a 
a ii aaa Ci ry) = 0 (9) 
(4) 


where: g = a,r, and B 


Three equations result from Lagrange’s equations which describe the motion of thy 
satellite. The three equations are: 


Cy ea ae 
tt (" cos’B e ) = 0 | {10) 
2 a (EL | a 
; de ( x) roosie( Ge) = -# 
d (2a 2 a 
7 ( oF )+? sin B cos p( -) = ap (12) 
From the equations of motion, two integrals result which are: 
22 dp ade 
r’ cos B OT constant (13) 
T+V = constant (14) 


Equation (13) results from integrating equation (10) and equation (14) simply states 
that the specific mechanical energy of the satellite remains constant. To change the in- 
dependent variable from ¢ to @, equations (1) , (2) , and (6) are used in conjunction 
with equation (13) and some initial conditions to form: 








dt _ r? cos i 
do- hy COS ig 


di 
1+ tan @-coti 16 ) (15) 


‘Letting w=p,/r, and using equation (15), the remaining equations of motion (11) 


- (12) can be rewritten as: 


2. 








di _ = 2Jusin@ cos @ sin i cos’ (16) 
dO ¢ cn 
+ 2Ju sin’@ cos*i 
cos i 
2 
I tu = \e cos*i— Je*u cos*i [2 on sin 6 cos @ (3 cos*i — 1) — u] 
— 4P te sin’@ cos & cos*i (3 — cos*i y [{c* + AJuc? sin? cost: (17) 
+ 4u*J? sin’@ cos*i } 
dQ _ = tan@ di 
a8 “sini dt (18) 
where: ¢ = COS% 
= SiN i 
3J,R? 
2p 


Equation (18) results from uncoupling the equations for Q and i. The angles Q 
and é can be uncoupled by applying the fact that the orbital p'ane must contain the 
velocity vector. The differential equations (16) - (18) are coupled by nonlinear terms and 
apparently cannot be solved analytically. If the right sides of equations (16) - (18) are 


expanded in a Taylor series expansion in powers of J, the equations simplify to: 
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Each of the negiccted terms in equations-(19) - (21) are indicated by the (O) symbols. 
These are terms which will be multiplied by / to the third power or higher. Equations 


(19) - (21) are identical to those used as a sta) ing point in the analysis of Eckstein, Shi 


and Kevorkian [ Ref. 6]. 


C. SOLUTION 


The initial breakthrough of an analytical solution of the equations of motion re- 
presented by equations (15) , (19) , (20) , and (21) was obtained by Danielson arid Snider 
[Ref. 1,7]. Further refinement of the solution was later formulated by Danielson, 
Sagovac, and Snider [Ref. 2, 8]. The thiee authors-developed the solution through the 
extensive use of an algebraic manipulation computer program called MACSY.1A. 
Through the use of an algorithm, MACSYMA was able to solve for the variables u ( 
orr),?,Q2,and dt/dé in terms of 6. Tc solution includes all terms multiplied ‘by J 
and excludes terms of order J? and higher. In 7. “er to maintzin a solution of order J 
when 0 < 1/J, the solutions also need to apprupriately include terms multiplied by J’0. 
The solutions which analytically demonstrate « relative accuracy of order_J’0 are: 
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D. SIMPLIFIED SOLUTION 


As shown in equations (22° - (73) , if @ is restricted such that @< 1/J, the solution 
should be of order J and the n2glected terms should be of order J*. For an Earth sat- 
ellite, J<3/2x 107, so for at least 100 revolutions the relative error should be on the 
order of 107°. If @ is restricted such that @<1, all of the terms of order J’@ in 
equations (22) - (26) can be neglected if a relative e:rer of order J is desired. By neg- 
lecting terms of order J’@, the solution simplifies considerably to: 
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E. THE CRITICAL INCLINATIONS 


As shown in equations (22) - (26) , the solution appears to be well bounded for al- 
most all inclinations. However, two particular inclinations immediately appear that may 
produce a singularity. A possible singularity occurs when the inclination is equal to ¢i- 
ther 63.43 degrees or 116.57 degrees. These two inclinations have produced mountains ; 
of literature and are well known as the critical inclinations. However, if the solution is 
replaced by the limit as the inclination approaches either critical inclination, the solution 
remains finite. More specifically, the solution at either critical inclination is as follows: 
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Clearly, equations (31) - (35) demonstrate that the solution is indeed finite for both 
critical inclinations. Equations (31)-- (35) are only valid for the critical inclinations and 
were first developed by Sagovac [Ref. 2] . The primary purpose in developing these 
equations was over a concern in computer programming. Some computers have major 
problems when a denominator approaches zero, and unlike humans, will not replace a 
solution with its limit. Therefore, depending upon the accuracy of a computer, 
equations (22) - (26) can replace equations (31) - (35) for inclinations near the critical 
inclinations. It should be noted, however, that the solution itself is valid and bounded 
for all inclinations. It is the limitation of the computer that creates the singularity. 

The simplified solution which is shown in equations (27) - (30) , is valid for all in- 
clinations. Since all terms of order J’@ have been neglected, the troublesome denomi- 
nators mentioned earlier do not appear. 


F. SPECIFIC MECHANICAL ENERGY 


For all satellites under the influence of conservative forces, the specific mechanical 
energy remains constant. Therefore, ar. ideal analytical check of the solution would be 
to sec if indeed the specific mechanical energy at any time is a constant. This simple 
check was performed by Danielson, Sagovac, and Snider [Ref 1, 2, 7, 8.] by substituting 
equations (22) - (26) into equations (7) and (8). The substitution yields: 


GM(1—e2) GAMJ,R7(1 — 3 sin”) 


aC 270 20 r(4) 1° 


+O, J’, ...) (36) 


The first two terms on the right side of equation (36) represent the initial specific 
mechanical energy. All other terms multiplied by J in equations (22) - (26) combine to 
zero when substituted into equations (7) and (8). Equation (36) demonstrates that by 
neglecting all terms of order J* and higher, the specific mechanical energy at any time 
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is precisely equal-to the initial specific mechanical energy. Obviously, the solution sat- 
isfies the requirement of constant specific mechanical energy. 











Ul. METHOD OF ANALYSIS 
A. ORBITAL PARAMETERS 


1. Argument of Latitude (6) 


Figure 2 illustrates that the position of a satellite at a particular time can be 
described by the argument of latitude (@) , the radius magnitude (r) , the inclination (i), 
and the longitude of the ascending node (Q). As shown in both the solution and the 
simplified solution, r, i, and Q are only functions of J and the argument of latitude (0) . 
Since J is a constant for all planets, a simple determination of these terms is trivial once 
6 is known. However, the determination of-6 is not trivial. Although it would be ideal 
for all of the equations to be analytical expressions, equations (26) , (30) , and (35) 
contain an integral that must be evaluated in order for @ to be determined. Herein lies 
the key to the solution. Given an elapsed time between observations, how can ¢ be 
precisely determined? Since the initial angular momentum (/,) is known, this term can 
be moved to the left side of equations (26) , (30) , and (35) to yield equations in the form 
of: 


0 
(t— tly = | (J, OL +f, 0)}a0 (37): 
‘Yo 


0 


If y was not a function of @, an evaluation of the right side of equation (37) 
could easily be conducted that would yield an analytic expression. However, r is also a 
function of @ and the only practical techni.ue in evaluating the integral is through nu- 
merical means. 

Several numerical methods could be used to evaluate the integral depending on 


the speed and accuracy one desires. Since accuracy and not speed is desired in this 


analysis, a Romberg integration routine was used to evaluate the integral. Since the 
right side of equations (26) , (30) , and (35) are sinusoidal in nature, the Romberg 
‘scheme converged quickly and accurately. 














Since @ defines the upper limit of the integral, in order to arrive at a solution, 
an initial @ must be estimated. Once @ is estimated, the integral can be numerically in- 
tegrated. and the result can be compared to the left side of equation (37). If the com- 
parison is accurate within some predetermined error, the iteration is complete and @ has 
been determined. If the comparison produces an error that is unacceptable, @ can be 
incremented either up or down and the integral can be reevaluated. Eventually, the it- 
eration will converge and @ will be determined. An algorithm of the iteration procedure 
is as follows: 


1. Estimate @. 

2. Evaluate the integral. 

3. Compare the result with the left side of equation (37). 
4. If outside the limit, go to (5). If within the limit, stop. 


5. Increment @ up or down as needed, go to (2). 


The determination of @ involves a combination of two errors. The first error is 
contained in the numerical evaluation of the integral itself, while the second error in- 
volves the comparison of the result of the integration with the left side of equation-(37). 
Unfortunately, the errors do not linearly combine, but rather multiply since the numer- 
ical evaluation of the integral is inherently nonlinear. In order to make the comparison 
error meaningful, the evaluation of the integral must be made as precise as possible. In 
order to avoid determining whether an error is due to computing or truncation errors, 
the numerical technique used in this analysis did not rely on a step size constraint. 
Therefore, the relative error, in general, can be specifically controlled. Since in this 
analysis, accuracy and not speed is desired, the Romberg integration technique was uti- 
lized. The Romberg technique does not depend on any specific step size and the evalu- 
ation of the integral is determined through a converging algorithm. Also, the relative 
error of the integration can be specifically controlled. In general, the relative error 
normally demanded in the integral evaluation was on the order of 107”, and the relative 
error of the comparison was on the order of 107”. Since the computer program utilized 
in the analysis was written for double precision accuracy, these types of relative errors 
presented no significant problems. The double precision accuracy enabled the computer 
program to calculate up to sixteen digit precision. 
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2. Radius Magnitude (7) 


From equations (22) , (27) , and (31) , it can be seen that the radius magnitude 
(r) is a function or J and 6. Oncé @ is known, r can be evaluated. From the appearance 
of equations (22) , (27) , and (31) , it is not obvious how r will behave as the orbit of a 
satellite progresses. However, from observations of actual satellite motion, it is clear 
that the orbit should behave elliptically with r varying from a minimum value at 
periapsis to a maximum value-at apoapsis. The magnitude of J plays an important role 
and fortunately for most planets, oblateness effects act as a perturbation in comparison 
to the main gravitational force. Therefore, a large value of J causes iarger variations. in 
r. Since equations (22) , (27) , and (31) contain a number of sine and cosine terms, a 
sinusoidal behavior should be expected. 


3. Inclination (7) 


The solution of the inclination is shown in equations (24) and (33) , and the 
inclination for the simplified solution is shown in equation (28). In general, these three 
solutions are quite similar. Again, once @ is known, i can be evaluated easily. It can 
be seen from equations (24) , (28) , and (33) , that / will vary slightly from an initial 
inclination as the orbit of a satellite progresses. Also, since a number of sine and cosine 
terms are present, the variation should be sinusoidal in nature. From inspection it is 
clear that the magnitude of the variation is dependent upon-the magnitude of J and the 
initial inclination. The variation of the inclination should not behave in a diverging 
fashion, but rather in an oscillatory fashion about some arbitrary mean inclination. This 
behavior is consistent with observations of actual measured satellite data. The driving 
factor in all inclination variations is the magnitude of J. Since for Earth, J, is on the 
order of 107°, these variations should be quite small. 


4, Longitude of the Ascending Node (2 ) 


The solution of the longitude of the ascending node (Q) is shown in equations 
(25) and (34) , and the longitude of the ascending node for the simplified solution is 
shown in equation (29). As expected, all solutions are quite similar. As with the case 
of r and i, Q can easily be determined once @ is known. Unlike the behavior of r and 
i, the variation of Q is very predictable and highly meaningful. With the presence of @ 
alone in equations (25) , (29) , and (34) , 2. possesses a linear relationship with 0 and 
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as @ increases with time, the variation of Q from Q, shouldbe linear. Depending upon 
the initial inclination, this variation will be either positive or negative. This type of be- 
havior is clearly consistent with the classical behavior known as nodal regression. For 
an oblate planet, nodal regression is a linear property whose magnitude and direction 
depends upon the radius magnitude and inclination of the-satellite. In equations (25) , 
(29) , and (34) , the radius magnitude is contained in the J term. Therefore, the mag- 
nitude of the nodal regression is entirely dependent upon the magnitude of J. From the 
analysis of the behavior of Q as @ increases, the nodal regression behavior should -be 
extremely obvious. 


B. ROMBERG INTEGRATION TECHNIQUE 


The Romberg integration technique is a powerful integration method in which ar- 
bitrary accuracy can be achieved in a relatively efficient manner. The method combines 
any type of relatively inaccurate quadrature method-with a Richardson extrapolation in 
order to quickly and accurately converge on a solution. In this analysis, a simple 
trapezoidal quadrature was initially used to estimate the integral and then a Richardson- 
extrapolation was used to improve the integration to the desired accuracy level. The 
trapezoidal quadrature first estimates the integral with a single interval. The estimate is 
then improved by using 2 intervals, 4 intervals, 8 intervals, etc. For purposes of iden- 
tification, the results can be labeled J} , J? , J; and so on. These results can be arranged 
in column form in preparation for a Richardson extrapolation and each new member 
represents the technique of halving the prior interval. The length of the column is de- 
termined by the accuracy that one desires. Once the first column is arranged, a 
Richardson extrapolation can be performed by the following equation. 


qiy” ss nf2 
t 4! 1 ( ) 


The values of J; can be arranged in tabular form as shown in Table 1. 














Table 1. A SCHEMATIC OF ROMBERG INTEGRATION 





To test for convergence, the value represented by J” is compared with the value 
represented by J”, . If these two values are within some predetermined error, then J?” 
becomes the evaluation of the integral. If convergence has not been reached, then an- 
other row is calculated and the process continues. An excellent example of the Romberg 
integration technique is shown in Ferziger [Ref.9]. In this example the following 
solution of the integral is desired. 


1 
I = | ea (39) 
0 


From elementary calculus, the exact solution is: 
Texact = 2-718281828 (40) 


The technique of Romberg integration of the integral is shown in Table 2. The rel- 
ative error of J} to If is 7.8110. The relative error of Jf to Ij is 1.97x 107". 
As can easily be seen, the integration is converging very nicely and the error found in the 


final solution is less than the error demanded within the Romberg integration scheme 
itself. 














Table 2. ROMBERG INTEGRATION 








re 1859140914 
1.753931092 | 


pa I. 727221905 
Fe 1. 720518592 1.718284155 1.718281842 1, 718281829 


The advantage of the Romberg integration technique over a simple quadrature 





method is obvious. The number of intervals that must be evaluated is very small and 
the relative accuracy is very high. In order to attain the accuracy that the Romberg 
technique delivers, the trapezoidal method would need to divide the integral into several 
more intervals. This would be highly inefficient. For smooth functions, the Romberg 
technique is very effective and efficient. Since equation (37) is sinusoidal in nature and 
thus relatively smooth, the Romberg integration technique was used to evaluate the in- 
tegral. If equation (37) had: not been so well behaved, another integration technique 
might-have been warranted. The Romberg integration scheme is the heart of the anal- 
ysis and can be found in the computer program shown in Appendix E. 














IV. METHOD OF COMPARISON 


A. NUMERICAL INTEGRATION COMPARISON 


In order to verify that the theory is valid for practical application, the solution must 
be compared-with proven numerical solutions and measured satellite data. By compar- 
ing the solution with a numerical integration of the equations of motion, theoretical ac- 
curacy can be specifically determined. As shown in equations (22) - (26) , theoretically 
the solution is accurate to order J’@. A numerical integration comparison will determine 
whether this prediction 1s correct. In order to verify the solution, the following param- 
eters will be compared. 


1. Delta-radius vector.(| Ar!) 
2. Earth arc angle (Y) 
3. Delta omega (AQ) 

Delta inclination (Af) 


4 
5. Delta theta (Aé@) 
6 
r | 





Delta radius relative error (| Ar| /r,) 
Delta theta relative error (A8/8,) 

8. Radial track error (RTE) 

9. Along track error (ATE) 

10. Cross track error (CTE) 


1. Delta Radius Vector 


A graphical representation of the delta radius vector (|Ar]) and the Earth arc 
length (WV) is shown in Figure 3. The delta radius vector is the magnitude of the vector 
separating the solution radius vector (r) from the numerical integration radius vector 
(r,). Mathematically, the delta radius vector can be expressed as: 


Ar = r-?, (41) 














The delta radius vector describes in overall terms-the global! error in the solution. 
Another common name for this error is the Euclidean normed difference in ephemerides. 
Although the delta radius vector provides ample information on the global error in the 
solution, this error can-also be expressed by a different parameter that will be called 
Earth arc angle. 








Figure 3. Delta radius vector and Earth arc angle. 


2. Earth Arc Angle 


Earth arc angle (‘V) is simply the angle between the two positions if viewed from 
sea level on Earth. For simplicity, the position at sea level was chosen such that the arc 
angle fram the center of the Earth was bisected. By using the law of sincs and cosines, 
the Earth arc angle is easily determined. Since satellites are tracked by instruments on 
the surface of the Earth, a better feel for the global error can be attained by determining 


the angle between the two positions. Most satellite tracking radars possess beamwidth 











and field of view limitations; therefore, Y will provide useful information on whether the 
solution is accurate enough fcr satellite tracking radars. 


3. Delta Omega, Delta Inclination, Delta Theta 


A. break down of the global error can be described in the errors of delta omega 
(AQ) , delta inclination (Ai) , and delta theta (A@) . As mentioned previously, AQ will 
provide an insight on the motion of the line of nodes, and specifically if nodal regression 
is present. The change in inclination. will provide information on the movement and 
stability of the orbital plane. The parameter in which all errors are based, A@ , will 
provide much information on the source of the global error. It is clear that small errors 
in A@ will contribute significantly to the accuracy of the solution. 


4. Relative Errors 


The verification of the solution will lie in the confirmation of the relative errors. 
The delta radius relative error and the delta theta relative error will demonstrate the ac- 
tual accuracy of the solution. Both parameters, Ar and A@ , demonstrate a theoretical 
error of J°@.. Therefore, the delta radius relative error should be on the order of J’0, 
while the delta theta relative error should be on the order of J’. Comparisons of the 
relative errors between the suiution and the numerical integration solution will provide 
the evidence for theoreticat confirmation. 


5. Track Errors 


Another me‘hod to break down the global error is in terms of track errors. 
Figure 4 shows a graphical representation of radial traci, error (RTE) , along track error 
(ATE) , and cross track error (CTE) . These errors can better be described by referring 
also to Figures 1 and 2. The radial track error is the error in the radial direction or in 
the B, direction. The along track error is the arc length error in the plane defined by the 
solution radius vector (r) and the B, direction. Together, these errors describe errors in 
three orthogonal directions or planes as compared with some reference position. The 
reference position in this case is the numerical integration solution. A mathematical 
derivation of these errors is as follows: 


Ar = rr, + Ar, 0 + A0, i+ AO, Q + AQ) — r(x, 8, i, Q) (42) 














Using equation (1), 
_ Ar = (r+ Ar(B, + AB) — 7B, 
Ar = r,B,+ArB,+ 1, AB, + ArAB, — 7,,B, 
Neglecting -higher order terms, 
Ar = ArB, +7,AB, 
Continuing, defining AB, , 
AB, = (B,-AB,)B, + (B,- 4B,)B, + (B3 + AB,)Bs 
Using the rotation transformation and after performing considerable algebra, 
B,-AB, = 0 


B,-AB, = (A0+AQ cos i) 


B;+AB, = (Aisin @ — AQ cos @ sin i) 
Therefore, using equation (45), 
Ar = (Ar)B, +7,(A0 + AQ cos JB, + r,(Ai sin 6 — AQ cos @ sin 1B; 
From equation (50), the track errors can easily be defined. 


RTE = r—-nr, 
ATE = »,(A@+ AQ cos i) 


CTE = 1,(Aisin @ — AQ cos @ sin i) 





(43) 


(44) 


(45) 


(46) 


(47) | 
(48) 


(49) 


(50) 











A graphical representation of the track errors is shown in Figure 4. 


ATE 





Figure 4. Track errors. 


Examination of equations (51) - (53) demonstrates that the track errors divide 
the global error into three distinct regions. Radial track errors obviously describe errors 
in the radial direction. Along track errors are similar in nature to Earth arc angle errors, 
but also include errors due to nodal regression. Cross track errors describe orbital plane 
errors in terms of both inclination errors and errors due to nodal regression. 

In general, all these parameters should give an excellent insight into the accu- 
racy of the solution. Also included in the numerical comparison will be the simplified 
solution and the two-body solution. The simplified solution has been previously pre- 
sented. The two-body solution can easilv be determined by simply setting J=0. The 
analysis of the numerical comparison wi!l demonstrate the strengths and weaknesses of 


all the different solutions. 
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B. MEASURED DATA COMPARISON 





In order to observe how well the solution models actual satellite motion, the sol- 
ution will be compared with actual measured data from operational satellites. To 
properly evaluate the solution, a wide range of orbit characteristics will be compared. 
These characteristics include orbits of various altitudes, inclinations, and eccentricities. 
Also included in this comparison will be the simplified solution, the two-body solution, 
and if available, some particular numerical solutions. 

The numerical solutions will consist of two forms. The first is a numerical solution 
that only includes perturbations involving J,,J;,J,, and J,. The second numerical 
solution will include the following perturbation effects. 


lL. J,,53,4,, and J, . 
. Atmospheric drag. 
. Sun gravitational effects. 


. Moon gravitational effects. 


wm Bb Ww Ww 


. Solar pressure effects. 


From the analysis, the accuracy of the solution and the simplified solution can be 
compared to a numerical solution as well as to actual measured data. The weakness of 
the two-body solution will also be demonstrated. In addition, the strengths of a well 
modeled numerical solution will clearly be seen. The identical error parameters described 
in the previous section will also be used in the measured data comparison. From this 
comparison, the advantages and disadvantages of the solution in regard to actual satcl- 
lite motion will clearly be demonstrated. 

















VY. RESULTS 


A. NUMERICAL INTEGRATION COMPARISON 


To verify that the theory is valid for practical applications, the solution was coin- 
pared with a proven numerical solution. A numerical integration computer program 
called UTOPIA is currently in use at -he Colorado Center for Astrodynamic Research 
(CCAR) located on the campus of the University of Colorado, Boulder, Colorado. 
UTOPIA is primarily used to model a wide range of perturbations and can predict sat- 
ellite motion with a high degree of accuracy. The UTOPIA computer program was de- 
veloped at the University of Texas, Austin, Texas, and is currently in use at several 
universities and research centers. The solution was compared with the UTOPIA sol- 
ution for a satellite with the following initial conditions: 


Yo = 7,386.18 kin 
ig = 90.03 degrees 
& = 0.003991 

@Qy = 224.38 degrees 
Q@ = 104.05 degrees 


Q, = 322.63 degrees 
hg = 54,205.18 kin’ ]s 
Po = 7,371.29 km 
Ig = 0.00 seconds 


In general, these initial cenditions represent a slightly retrograde orbit of small ec- 
centricity at an altitude of approximately 1000 kilometers. Essentially, it is a polar orbit 
at an altitude where several satellites are currently in motion. From the initial condi- 
tions, the orbit should demonstrate a slight easterly nodal regression. But, since the in- 
clination is so close to 90 degrees, some integration routines might predict zero nodut 
regression. In this comparison, UTOPIA only modeled the J, perturbation; therefore, 
the solution should compare well if the theory is valid. All error parameters depicted in 
the comparison were calculated in the following manner. 


A(ErrorParameter) = Theoretical Solution — UTOPIA Numerical Solution (34) 
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All solutions were compared at one hour intervals over two separate periods of time. 
One comparison is for a time period of one day while the second is for a time period of 
three days. The three day comparison was constructed to illustrate the effect of long 
term errors while the one day comparison allows for a more detailed analysis during the 
first few hours of motion. The period of rotation for the satellite is about 105.26 minutes 
which equates to approximately 13.68 orbits per day. The results of the numerical sol- 


ution comparison are shown in Appendix A. 


1. Delta Radius Vector Comparison 


The comparison of the delta radius vector is shown in Figures 5, 6, and 7 in 
Appendix A. Figure 5 includes the comparison of the soJution, the simplified solution, 
and the two-body solution to the numerical solution. If the solution matched the nu- 
merical solution exactly, the delta radius vector would be zero. As shown in Figure 5, 
the solution and the simplified solution match extremely well with the numerical solution 
while the two-body solution contains gross errors. 

A more detailed plot of the delta radius vector comparison is shown in Figures 
6 and 7. In these figures, the two-body solution is excluded. The difference in the sol- 
ution and the simplified solution can clearly be seen. The simplified solution produces 
a diverging sinusoidal response about the solution. However, up to approximately four 
hours of motion, the solution and the simplified solution are nearly identical. The 
sinusoidal behavior of the simplified solution can be attributed to the fact that all terms 
multiplied by J’°0 have been neglected. As 6 grows with time, these terms become sig- 
nificant in the solution. As shown specifically in Figure 7, the average delta radius vec- 
tor of the simplified solution clearly diverges from the solution. 

Figures 5, 6, and'7 also demonstrate that for at least one day, the delta radius 
vector for the solution and the two-body solution are nearly linear as a function of time. 
Other comparisons will determine whether this relationship holds true and will be shown 
later. As mentioned earlier, the delta radius vector is a global error. As shown in Fig- 


ures 5, 6, and 7, the solution compares well globally with the numerical solution and 


demonstrates a great improvement over the two-body solution. 














2. Earth Arc Angle Comparison 


The comparison of the Earth arc angle is shown in Figures 8, 9, and 10 in Ap- 
pendix A. From inspection, these plots are nearly identical in appearance to the delta 
radius vector plots. This is expected since both the delta radius vector and the Earth are 
angle represent global errors. Figure 8 clearly illustrates the large error generated by the 
two-body solution. The two-body solution produces unsatisfactory. long term satellite 
position prediction. After one day, a tracking radar would have a difficult time detecting 
a satellite with a position error of over 80 degrees. 

Figures 9 and 10 present much more encouraging results. Again, the simplified 
solution responds in a sinusoidal behavior about the solution. After one day, the posi- 
tion error of the solution is only approximately 0.15 degrees. Clearly, the solution and 
the simplified solution are superior to the two-body solution. Most tracking radars can 
easily handle daily position errors of 0.15 degrees. In general, the solution and the sim- 
plified solution agree very well with the numerical solution. 


3. Delta Omega Comparison 


The comparison of the delta omega angle is shown in Figures 11 and 12 in Ap- 
pendix A. At first glance, the solution and the simplified solution in Figure 11 appear 
not to agree well with the numerical solution. However, the scale of delta omega is 
multiplied by 10™. The numerical solution parallels the two-body solution nearly ex- 
actly and predicts almost no change in 22. In other words, the numerical solution pre- 
dicts no nodal regression. Easterly nodal regression is represented by a positive delta 
omega; therefore, it is clear that the solution and the simplified solution predict greater 
nodal regression than the numerical solution. 

On a larger scale, all three solutions are essentially identical. Since the initial 
inclination is so close to 90 degrees, small discrepancies are not surprising. The delta 
omega plot does, however, invoke confidence in the solution. Although the initial in- 
clination is very close to 90 degrees, the solution, the simplified solution, and the nu- 
merical solution predict easterly nodal regression. This result is significant. Even initial 
inclinations close to 90 degrees produce nodal regression in the correct direction for the 


solution, the simplified solution, and the numerical solution. 








4, Delta Inclination Comparison 


The comparison of the delta inclination angle is shown in Figures 13 and 14 in 
Appendix A. On a larger scale, all of the solutions compare very well. On the scale 
shown in Figure 10, the two-body solution and the numerical solution are nearly iden- 
tical. The solution and the simplified solution oscillate about an error of approximately 
—2.0 x 107° degrees. Obviously, this error is extremely small. In general, the solution 
and the numerical solution agree very well. 

An interesting aspect of the delta inclination comparison is the sinusoidal be- 
havior of the solution and the simplified solution. This type of behavior is precisely what 
was predicted in the earlier analysis. Figure 14 demonstrates that this behavior contin- 
ues for even longer periods of time. On a larger scale, this type of motion would not 
be detectable. 


5. Delta Theta Comparison 


The comparison of the delta theta angle is shown in Figures 15, 16, and 17 in 
Appendix A. This comparison confirms the results found in the earlier comparisons. 
The two-body solution produces very large errors, while the solution and the simplified 
solution agree very well with the numerical solution. Figures 16 and 17 again illustrate 
the typical sinusoidal response of the simplified solution about the solution. Since the 
delta theta error produces all other errors, the excellent results found in the earlier 
comparisons are now not surprising. 


6. Delta Theta Relative Error Comparison 


The comparison of the delta theta relative error is shown in Figures 18, 19, and 
20 in Appendix A. As shown in Figure 18, the two-body solution demonstrates a rela- 
tive error of 2.3/, while the relative errors of the solution and the simplified solution are 
much smaller. In more detail, Figures 19 and 20 indicate that the relative error of the 
solution is 2.8/?. This result confirms the theoretical prediction that .ne delta theta 
relative error of the solution would be on the order of J’. Figures 19 and 20 illustrate 
that initially the relative error of the simplified solution is also on the order of J*. But 
as @ increases with time, the relative error grows in a sinusoidal fashion. This result is 


expected since the simplified solution neglects all the terms multiplied by J’0. In gen- 


3] 











eral, the results shown in Figures 18, 19, and 20 confirm the theoretical relative accuracy 
of delta theta that was predicted-in the earlier analysis. 


7. Delta Radius Relative Error Comparison 


The comparison of the delta ratlius relative error is shown in-Figures 21, 22, and 
23 in Appendix A. As shown in Figuré 21, the two-body solution produces a relative 
error that is“linear in time-and proportional to 2.3/8. Again, the relative errors of the 
solution and the simplified-solution are magnitudes smaller. 

Figures 22 and 23 present in more detail the relative errors of the solution and 
the simplified solution. The relative error of the solution is very linear and proportional 
to 2.80. The relative error of the simplified solution is sinusoidal in nature and di- 
verges from the solution. However, for up to four hours of motion, the relative error 
of the solution and the simplified solution are nearly indistinguishable. Again, the re- 
sults from this comparison confirm the theoretical prediction that the delta radius rela- 
tive error of the solution would be on the order of J%. 


8. Radial Track Error Comparison 


The comparison of the radial track error is shown in Figures 24, 25, and 26 in 
Appendix A. As shown in Figure 24, the two-body solution oscillates about an error 
of approximately —11.0 kilometers, while the solution and the simplified solution both 
produce errors that are dramatically. smaller. From inspection, the two-body solution 
also appears to be slowly converging as time increases. 

In Figures 25 and 26, the solution and the simplified solution produce con- 
trasting behaviors. While the solution remains relatively constant, the simplified sol- 
ution slowly diverges from zero. These two different responses continue even after three 
days of motion. Again, the neglected J’@ terms cause the significant divergence of the 
simplified solution. Not surprisingly, the solution and the simplified solution are clearly 
superior to the two-body solution. : 


9, Along Track Error Comparison 


The comparison of the along track error is shown in Figures 27, 28, and 29 in 
Appendix A. The results presented in this comparison parallel the results found in the 
earlier comparisons. Since the inclination is so close to 90 degrees, the AQ contribution 
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is negligible and the A@ contribution strongly influences the responses. As a result, the 
along track error comparison is practically a mirror image of the delta theta comparison. 


10. Cross Track Error Comparison 


The comparison of the cross track error is shown in Figures 30 and 31 in Ap- 
pendix A. The cross track error is strongly influenced by Ai and AQ. Since the two- 
body solution produced good results with these two parameters, it is not surprising that 
the two-body solution agrees well with the numerical solution. Fortunately, the errors 
produced by the solution and the simplified solution are also very small. The solution 
produces a maximum cross track-error of approximately + 0.5 kilometers after one day 
of motion, and approximately + 1.3 kilometers after three days of motion. 

Clearly in this comparison, the two-body solution is superior. However, the 
large errors produced by the two-body solution in the other comparisons easily over- 
whelm these results. In global terms, the two-body solution is no match for either the 
solution or the simplified solution. 


B. MEASURED DATA COMPARISON 


The solution -was compared with actual measured satellite data to determine the al- 
titude band where the theory works best. The measured satellite data was obtained from 
the First Satellite Control Squadron (1SCS) located at Falcon Air Force Base, Colorado. 
The First Satellite Control Squadron tracks several satellites for the Air Force and was 
able to supply measured data for three separate satellites. The three satellites are cur- 
rently in motion and occupy orbits that are labeled Near Earth, Semisynchronous, and 
Geosynchronous, respectively. All error parameters compared in the earlier numerical 
comparison were also compared in this comparison using the measured data as a refer- 
ence. 

Included in all the comparisons were the solution, the simplified solution, the two- 
body solution, and two numerical solutions. The two numerical solutions were also 
supplied by the First Satellite Control Squadron and are labeled Spacom | and Spacom 
2, respectively. The Spacom | solution includes all perturbation effects, while the 
Spacom 2 solution only includes the Earth’s harmonic perturbations. All error param- 
eters in this comparison were calculated in the following manner. 
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A(ErrorParameter) = Test Solution — Measured Data Solution (55) 


Unfortunately, the First Satellite Control Squadron only records measured data 
when an update of their numerical solution is required. Routine updates are usually 
conducted after about seven days of motion. Therefore, satellite data for one month 
usually consists of only four data points. Although-more data points are needed for a 
more detailed analysis, a long term analysis can still be conducted. The analysis of each 
type of orbit will be presented separately. 


1. Near Earth Orbit Comparison 


The near Earth orbit comparison possesses the following initial conditions. 


Y = 7,776.58 km 

in = 98.81 degrees 

€y> = 0.0003071 

@) = 9.57 degrees 

0) = 149.14 degrees 

Qy9 = 37.10 degrees 

lg = 53,664.37 kin?|s 

Po = 7,224.89-km 

fy = 0000Z 26 July 1990 


The initial conditions of this satellite represent a retrograde orbit of small ec- 
centricity at an altitude of approxiinately 850 kilometers. The period of rotation for the 
satellite is about 101.89 minutes which equates to approximately 14.13 orbits per day. 
From the initial conditions, J, should be the dominant perturbation. The orbit should 
demonstrate noticeable easterly nodal regression. If the theory is valid, both the sol- 
ution and the simplified solution should agree well with the numerical solutions and the 
measured data. 

The results of the near Earth orbit comparison are shown in Figures 32 - 43 in 
Appendix B. As shown in the figures. the solution and the simplified solution agree very 
well with both the Spacom | solution and the measured data. The fact that the solution 
and the simplified solution produce such excellent results verifies that J, is the dominant 
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perturbation: ‘for this satellite. Figures 32 - 43 also demonstrate the larger errors 
produced by the two-body solution. In almost every comparison, both the solution and 
the simplified solution are far superior to the two-body solution. 

One surprising result is the poor comparison produced by the Spacom 2 sol- 
ution. In every comparison the Spacom 2 solution either models the two-body solution 
exactly or produces results that are inferior to the two-body solution. It is clear that the 
Spacorn 2 solution doesnot model the Earth’s harmonic forces correctly. An explana- 
tion for the poor results cannot be determined in this analysis. A detailed analysis of the 
force modeling in the Spacom 2 solution must be completed in order to adequately ex- 

- plain the unsatisfactory results. 

The delta omega comparison in Figure 34 demonstrates the easterly nodal re- 
gression produced by the solution, the simplified solution, the Spacom 1 solution, and 
the measured data. The two-body solution represents zero nodal regression. Figure 35 
presents the-delta omega comparison at a much smaller scale and excludes the two-body 
solution. In this figure, much more detail can be observed. 


There is only one comparison in which the results are mixed. The radial track 
error comparison in Figure 40 indicates that the solution produces a small improvement 
over the two-body solution while the simplified solution actually produces a greater er- 
ror. In comparison with the along track errors, these errors are small. It is interesting, 
however, that the radial track error comparison produces such mixed results. In general, 
both the solution and the simplified solution produce results that are in excellent agree- 
ment with the measured data for this near Earth satellite. 














2. Semisynchronous Orbit Comparison 


The semisynchronous orbit comparison possesses the following initial condi- 
tions. 


’ = 26,407.70 km 
in = 63.66 degrees 
€y = 0.005860 

Wy = 318.19 degrees 
Q) = 328.49 degrees 
Qo = 92.13 degrees 


hg = 102,892.59 km?/s 
Po = 26,559.96 kin 
fy = 0000Z 22 March 1990 


The initial conditions of this satellite represent a direct orbit of small eccentricity 
at an altitude of approximately 20,000 kilometers. The period of rotation for the satellite 
is about 717.96 minutes which equates to approximately 2.01 orbits per day. An im- 
portant aspect-of the orbit is that-the initial inclination is-very close to the critical incli- 
nation of 63.43 degrees. Although the initial inclination is not exactly that of the critical 
inclination, an evaluation of the solution and the simplified solution near this important 
inclination can be made. From the initial conditions, the orbit should demonstrate 
substantial westerly nodal regression. Also, at this altitude, the dominance of the J, 
perturbation should be diminished. Other perturbations that are not modeled should 
make a considerable contribution to the errors in the comparison. If the theory is valid, 
both the solution and the simplified solution should show a great improvement over the 
two-body solution. 

The results of the semisynchronous orbit comparison are shown in Figures 44 - 
53 in Appendix C. As predicted earlier, the solution and simplified solution produce 
results that are superior to the results produced by the two-body solution. Figures 44 
and 45 present the global errors of all the solutions. In global terms, the solution and 
the simplified solution reduce the error of the two-body solution by nearly one half. In 
effect, the J, perturbation accounts for approximately one half the error produced by the 
two-body solution. The remaining error which is represented by the solution and the 
simplified solution is caused by other perturbing forces. Unfortunately, the results of the 
Spacom 2 solution were not available. 
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The delta omega comparison in Figure 46 demonstrates the easterly nodal re- 
gression produced by the solution, the simplified solution, the Spacom 1 solution, and 
the measured data. Again, the two-body solution represents zero nodal regression. It 
is clear that at this altitude, the J, perturbation produces the majority of the nodal re- 
gression. The delta inclination comparison in Figure 47 indicates that the solution and 
the simplified solution produce results that are not much better than the results 
produced by the two-body solution. However, the error after 30 days of motion is ex- 
tremely small. On a larger scale, the solutions would seem identical. Since the inclina- 
tion is very near the critical inclination, these results produce more evidence in support 
of the theory. Clearly, the solution and-the simplified- solution are bounded at this in- 
clination. The delta.theta comparison in Figure 48 demonstrates that the majority of the 
error produced by the solution and the simplified solution originates in the delta theta 
error. It is clear that the two-body solution underestimates the value of 6 while the 
solution and the simplified solution overestimate the value of 6 . 

The relative error comparisons are shown in Figures 49 and 50. While the delta 
theta relative error for the solution, the simplified solution, and the two-body solution 
is approximately 15.0 x 107%, the relative error produced by the Spacom | solution is far 
superior. This result is expected since the Spacom 1 solution models several more in- 
fluential perturbations. The delta radius relative error comparison again demonstrates 
in global terms the amount of improvement that the solution and the simplified solution 
provide over that of the two-body solution. 

The track error comparisons in Figures 51 - 53 produce mixed results. While 
the two-body solution produces less radial and along track errors, the solution and the 
simplified solution produce much less cross track error. In comparison with the along 
track and cross track errors, the radial track errors are small. The poor results produced 
by the solution and the simplified solution in the-along track error comparison is due 
primarily to the large error in AO. The very large error produced by the two-body sol- 
ution in the cross track error comparison is due primarily to the very large error in AQ. 

In summary, although the solution and. the simplified solution are superior to 
the two-body solution, the Spacom | solution models the satellite motion more precisely. 
However, the primary reason that the solution and the simplified solution are superior 
to the two-body solution is due exclusively to a better modeling of nodal regression or 
the angle 2. It is clear that the solution and the simplified solution model the J, per- 
turbation extremely well. The Spacom I solution is expected to perform better since it 


models more perturbing forces. 

















3. Geosynchronous Orbit Comparison 


The geosynchronous orbit comparison possesses the following initial conditions. 


% = 42,156.57 km- 

io = 1.09 degrees 

€& = 0.0002341 

Wy = 320.06 degrees 

6) = 331.32 degrees 

Qy = 334.85 degrees 

No = 129,644.14 km/s 
Po = 42,166.25 km 

fy = 0000Z 21 July 1990 


The initial conditions of this satellite represent a direct orbit of small eccentricity 
at an altitude of approximately 35,800 kilometers. The period_of rotation for the satellite 
is about 1436.69 minutes which-equates to approximately 1.00 orbit per day. Since the 
‘initial inclination ‘is slightly greater than zero, the orbit should demonstrate westerly 
nodal regression. However, since the altitude is so large, other perturbing forces that 
are not modeled may influence nodal regression. At a geosynchronous altitude, the 
magnitude:of other perturbing forces approach that of J,. Since at this altitude the ef- 
fect of J, is so diminished, some comparisons of the solution, the simplified solution, and 
the two-body solution may be nearly identical. As a result, the theory may not be any 
better than the two-body theory for satellites in a geosynchronous orbit. 

The results of the geosynchronous orbit comparison are shown in Figures 54 - 
63 in Appendix D. The global error comparisons are shown Figures 54 and 55. In 
global terms, the solution and the simplified solution produce results that are surpris- 
ingly superior to the results produced by the two-body solution. Evidently, for this sat- 
ellite, the J, perturbation jis still quite dominant. However, the other comparisons may 
present a different picture. Once again, the Spacom 2 solution generates very poor re- 
sults, 

The delta omega comparison in Figure 56 indicates that the actua! perturbing 
forces produce easterly nodal regression. Conversely, the solution and the simplified 
solution predict westerly nodal regression. It is obvious that other perturbing forces 


38 














influence the nodal regression of this satellite. Although the Spacom | solution is su- 
perior, even this accurate numerical solution has trouble predicting the value of Q. The 
solution and the simplified solution also produce poor results in the delta inclination 
comparison in Figure 57. All solutions, except for the Spacom 1 solution, produce 
identical results. Again, on a larger scale, all of the solutions would seem nearly identi- 
cal. However, this detailed analysis does demonstrate a weakness in the theory. The 
delta theta comparison in Figure 48 indicates that the solution and the simplified sol- 
ution are inferior to all solutions including the two-body solution. Clearly, other per- 
turbing forces are at work. 

The relative error comparisons are shown in Figures 59 and 60. The delta theta 
relative error comparison simply repeats the results found in the delta theta comparison. 
However, the delta radius relative error comparison is much more reassuring. Again, in 
global terms, the solution and the simplified solution produce better results than the 
two-body solution. 

The track error comparisons in Figures 61 - 63 produce mixed results. The ra- 
dial track error comparison indicates that initially the Spacom 1 solution is inferior to 
all other solutions. However, after 21 days of motion, Spacom | is the superior solution. 
Once again, the radial track errors are small when compared to the along and cross track 
errors. The clue to the favorable global results of the solution and the simplified solution 
is found in the along and cross track error comparisons. The solution and the simplified 
solution perform much better than the two-body solution in the along track error com- 
parison. Although the two-body solution is superior to the solution and the simplified 
solution in the cross track error comparison, the difference is small. It is clear that the 
solution and the simplified solution are superior to the two-body solution due to a much 
smaller along track error. 

In summary, although the solution and the simplified solution are superior to 
the two-body solution, other perturbing forces greatly influence the satellite’s motion. 
At this altitude, the solution and the simplified. solution simply do not model the satel- 


lite’s motion well. Other perturbing forces must be modeled at this altitude if proper 
satellite position prediction is desired. 

















VI. CONCLUSIONS.AND RECOMMENDATIONS 


An analysis was-conducted on a perturbation solution of the main problem in arti- 
ficial satellite theory. The purpose of the analysis was to compare the solution with 
proven numerical solutions and actual measured satellite data in order to determine if 
the theoretical work is valid and practical. From the analysis, the following conclusions 
can be made. 


1. The solution and the simplified solution are both significantly more accurate than 
the two-body solution. The relative error of the two-body solution is on the order 
of J@é While the relative error of the solution and the simplified solution is on the 
order of J‘6. 


2. The real physical effects of the orbit are easily distinguishable in both the solution 
and the simplified solution. 


3, The solution and the simplified solution compare extremely well with a proven 
numerical solution for at least 41 revolutions with a relative error on the order of 
J. 


4, The solution and the simplified solution compare extremely well with actual meas- 
ured satellite data for at least 297 revolutions at altitudes where the J, perturbation 
dominates ( e.g., near Earth orbits ). For a satellite in orbit at an altitude of around 
1000 kilometers, the solution and the simplified solution reduce the error of the 
two-body solution by approximately 95%. 


5. The solution and the simplified solution compare less favorably with actual meas- 
ured. satellite data at semisynchronous and geosynchronous altitudes. At these al- 
titudes, however, the solution and the simplified solution reduce the error of the 
two-body solution by at least 50%. 


6. The solution and the simphified solution are free of singularities and are valid for 
all orbital parameters. 


Clearly, the solution and the simplified solution model the J, perturbation very well. 
The equations are easy to implement and can provide quick and accurate predictions of 
satellite motion. However, other types of analytical solutions exist that are more accu- 
rate than the solutions described here. 

One such solution was developed by Coffey and Alfriend [Ref. 10] through re- 
search that was conducted by Dep.:. CRef. 11], Coffey and Deprit [Ref. 12], and 
Alfriend and Coffey LRef. 13] . The solution is called the Analytic Orbit Prediction 
Program gencrator or (AOPP). Although the program is very accurate, AOPP exten- 


40 











sively utilizes four different Hamiltonian transformations. As a result, the real physical 
effects of the orbit are not easily distinguishable. 

The beauty of the solution and the simplified solution is their similarity in form to 
the well known two-body solution and the fact that a satellite’s position can be easily 
predicted by evaluating only one integral. Once @ has been determined, all other orbital 
parameters can be calculated easily. The structure of the solution and the simplified 
solution is ideal for implementation with onboard spacecraft computers. 

Before the solutions can be adapted for practical applications, more examination 
and testing of the theory is required. In order to provide more confidence in the theory, 
the following recommendations are suggested. 


1. The solution and the simplified solution need to be compared-to a numerical inte- 
gration of the equations of motion for at least 100 revolutions to confirm the the- 
oretical accuracy for long term satellite motion. 


2, The solution and the simplified solution need to be compared to several more di- 
verse sets of actual measured satellite data. 


3. To increase precision, the solution needs to include the higher order zonal har- 
monics of the gravitational potential ( e.g., J; , Js , Js, ete. ). 


4. For spacecraft computer implementation, the Lagrangian coefficients of the state 
transition matrix need to be determined. 


For onboard spacecraft navigation, computers make use of the state transition ma- 
trix. Currently the Lagrangian coefficients of the two-body solution are the only matrix 
elements that have been determined. An excellent formulation of the two-body state 
transition matrix is shown by Battin [Ref. 14]. Once the Lagrangian coefficients of the 
solution are developed, onboard spacecraft navigation can be greatly improved. 
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PROGRAM COL02 


RrekeketetedeKA ATK IKI RI 


ve * 
* MAIN PROGRAM % 
* * 


kikkkkkkkkk kiki 


IMPLICIT DOUBLE PRECISION (A-I,M-Z) 
CHARACTER*20 LINE 
DIMENSION M( 100) ,MD(100) ,E( *00) ,WC 100) ,WD( 100) ,OM( 100) ,OMD( 100) 
DIMENSION I( 100) ,ID(100) ,F(1v1) ,FD( 100), EC(100), ECD( 100), A( 100) 
DIMENSION R(100),H( 100), N(100), TH( 100) , THD( 100) 
DIMENSION RF(100), 1F(100), IFD(100), OMF( 100) , OMFD( 100) , THF( 100) 
DIMENSION THED( 100) , P(100), JORBIT( 100), DR(100), DID(100), DTHD(.100) 
DIMENSION DOMD(100), RX( 100), RY( 100) ,RZ(100), RFX(100), RFY( 100) 
DIMENSION RFZ(100), DRV( 100), ARG( 100), ARCD( 100), DAY( 100), HX(100) 
DIMENSION HY(100), HZ(100) , VxX(1003 ,V¥(100) , vZ(100), DT(100), NX( 100) 
DIMENSION NY( 100) ,N2(100), RDV( 100), EX(100), EY(100), E2( 100) 
DIMENSION NDE(100), EDR(100), V(1003,HT(100), RDRF( 100) , INTA( 100) 
DIMENSION ROMA( 100) , TH. (100), ATE( 100), CTE( 100) 
COMMON/OBLATE1/DAY ,RX_RY,RZ, vx, VY,VZ,DT,HX,HY ,HZ,NX Os 
abi pe eae R,¥ SEX, FL, EZ, MU, NDE, EDR, H, N, E,P,1,0M,W,F 
COMMON/OBLATE3/PI,EC,M,4 JHE,ER, TH, THD, RTD, MD, WD, OmD ,1D, ECD 
cern eta aaa LINE, J,THE.,THFD,IF, TFD ,OMF, OMED ,RF, INT, ROM 
COMMON/OBLATES/RJ,DR, DID, DTHD, DOMD, ESTERR, ACTERR, TERROR, JVER 
COMMON/OBLATE6/RFX, RFY, REZ, ARC, ARCD, RDRF, DRV, RJ2,5N, JORBIT 
COMMON/OBLATE7/INTA, ROMA, THO, ATE, CTE 


PRINT*,'ENTER VERSION ( SOLUTION = 1, SIMPLE = 2, TWO BODY 
READ* , JVER 
IF(JVER. EQ. 1. OR. JVER. EQ. 2. OR. JVER. EQ. 3) THEN 
GOTO 20 
ELSE 
GOTO 10 
ENDIF 


" 
LS’) 
~~ 


PRINT*,'ENTER FIRST POINT’ 
READ* ,K 

PRINT*,'ENTER FINAL POINT’ 
READ* , KK 


PRINT*, ENTER RJ2' 
READ" , RJ2 
IF(RJ2. EQ. 1. ODO) THEN 

RJ2=0. 00108263D0 
ENDIF 


LINE=' -------------------- ; 
PI=3. 141592653589793238462643D0 
RTD=180. OD0/PI 

ER=6378. 137D0 

HTS=1. 0D0/806. 812D0 

MU=3. 98600436D5 


OPEN(UNIT=2, STATUS='OLD', FILE='/COLO2 OUT A') 
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OPEN(UNIT=3, STATUS='OLD', FILE='/COLO2 DSS A‘) 
OPENCUNIT=4, STATUS='0OLD', FILE='/COLO20 DSS B') 


CALL DATA 
CALL ELEMENTS 
CALL PINITIAL 


J=1 . 
WRITE(6,6000) ‘POINT = ',J 
WRITE(6,6000) ‘INTEGRATE COMPLETED' 
WRITE(6,6000). ‘FORMULA COMPLETED’ 
WRITE(6,6000}) 'INERTIAL COMPLETED' 


DO 30 J=K, KK 
WRITE(6,6000) ‘POINT = ',J 
CALL INTEGRATE 
WRITE(6,6000) ‘INTEGRATE COMPLETED' 
CALL FORMULA 
WRITE(6,6000) ‘FORMULA COMPLETED' 
CALL INERTIAL 
WRITE(6,6000) ‘INERTIAL COMPLETED' 
CONTINUE 
CALL RESULTS 


CLOSE( UNIT=2) 
CLOSE(UNIT=3) 
CLOSE( UNIT=4) 


FORMAT(3X,A,13) 
STOP 
END 


dedekdedetesededeedetetekieiekdevenietesitesede 
ate ’, 
fe ve 
% SUBROUTINE DATA % 
we 
vs 


detededeteacdievedesieteseleleteek kes sess 


SUBROUTINE DATA 

IMPLICIT DOUBLE PRECISION (A-I,M-Z) 

CHARACTER*20 LINE 

DIMENSION N(100) ,MD( 100) ,E( 100) ,W( 100) ,WD( 100) ,OM( 100) , OND( 100) 
DIMENSION I( 100) ,ID(100) ,F( 100) ,FD( 100) ,EC( 100) ,ECD( 100) ,A( 100) 
DIMENSION R(100) ,H( 100) ,N( 100) , TH( 100) , THD( 100) 

DIMENSION RF( 100) ,IF( 100) ,IFD( 100) ,OMF( 100) ,OMFD(100) , THF( 100) 
DIMENSION THFD(100) ,P( 100) , JORBIT( 100) ,DR(.100) ,DID( 100) , DTHD( 100) 
DIMENSION DOMD( 100) ,RX( 100) ,RY( 100) ,RZ( 100) ,RFX( 100) ,RFY( 100) 
DIMENSION RFZ( 100) ,DRV( 100) ,ARC( 100) ,ARCD( 100) , DAY( 100) , HX( 100) 
DINENSION HY(100) ,H2( 100) , VX( 169) ,V¥( 100) ,VZ( 100) ,DT( 100) ,NX(100) 
DIMENSION NY( 100) ,NZ( 100) , RDV( 100) ,EX(100) , EY( 100) ,EZ( 100) 
DIMENSION NDE( 100) ,EDR(100) ,V( 100) ,HT( 100) ,RDRF( 100) , INTA( 100) 
DIMENSION MONTH( 100) ,DATE( 100) , HOUR( 100) ,MIN( 100) ,SEC( 100) 
DIMENSION ROMA( 100) , THJO( 100) ,ATE( 100) ,CTE( 100) 
CONMON/OBLATE1/DAY ,RX,RY,RZ,VX,VY,VZ,DT,HX,HY,HZ,NX,NY,NZ,K,KK 
CONMON/OBLATE2/RDV,R,V,EX,EY,EZ,MU,NDE,EDR,H,N,E,P,1,OM,W,F 
COMMON/OBLATE3/FPI,EC,M,A,HT,ER,TH,THD,RTD,MD,WD,OMD,ID,ECD 
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agagggaNqa 





COMMON/OBLATE4/FD, LINE ,J,THF ,THFD, IF, IFD,OMF,OMFD,RF, INT ,ROM 
COMMON/OBLATES/RJ,DR,DID,DTHD,DOMD ,ESTERR, ACTERR, TERROR, JVER 
COMMON/OBLATE6/RFX,RFY,RFZ,ARC,ARCD,RDRF ,DRV,RJ2,JN,JORBIT 
COMMON/OBLATE7 /INTA, ROMA, THJO, ATE , CTE 


READ IN EMPHERIS DATA 
OPEN(UNIT=1, STATUS='OLD’, FILE='/COLO2 DAT’, ACTION='READ') 


DO 10 J=1, KK 
READ(1,*) MONTH(J) ,DATE(J) ,HOUR(J) ,MIN(J) ,SEC(J) 
READ(1,*) RX(J),RY(J) ,RZ(J) 
RX(J)=RX(J) /1000. ODO 
RY(J)=RY(J)/1000. ODO 
RZCJ)=RZ( J) /1000. ODO 
READ(1,*) VX(J),VY(J),VZ(J) 
VX( J)=VX(J) /1000. ODO 
VY(J)=VY(J)/1000. ODO 
VZ(J)=VZ(J)/1000. ODO 

CONTINUE 


CONVERT PARAMETERS 


DO 20 J =1, KK 
DAY( J)=DATE( J)+( (3600. ODO*HOUR( J)+ 
(60. ODO*MIN( J)+SEC(J))))/86400. ODO 
DT( J)=(DAY( J) -DAY( 1) )*24. OD0*3600. ODO 
CONTINUE 


CLOSE(UNIT=1) 


RETURN 
END 


dededesvietetekdeksverhkseekKkK IRI ase 


ste Se 
rh rh 


* SUBROUTINE ELEMENTS * 
* % 


dededksedevereredesksovevedesetereiesesedeteveesee ees 


SUBROUTINE ELEMENTS 

IMPLICIT DOUBLE PRECISION (A-I,N~-Z) 

CHARACTER*20 LINE 

DIMENSION M(100) ,MD( 100) ,E( 100) ,WC 100) ,WD( 100) ,OM( 100) ,OMD( 100) 
DIMENSION I(100) ,1D( 100) ,F(100) ,FD( 100) ,EC( 100) ,ECD( 100) ,A( 100) 
DIMENSION R(100) ,H( 100) ,N( 100) ,TH( 100) , THD( 100) 

DIMENSION RF( 100) , IF( 100), IFD(100) ,OMF( 100) ,OMFD( 100) , THF( 100) 
DIMENSION THFD( 100) ,P( 100) , JORBIT( 100) ,DR( 100) ,DID( 100) , DTHD( 100) 
DIMENSION DOMD(100) ,RX( 100) ,RY( 100) ,RZ( 100) ,RFX(100) ,RFY( 100) 
DIMENSION RFZ( 100) ,DRY( 100) ,ARC(100) , ARCD( 100) , DAY( 100) , HX( 100) 
DIMENSION HY( 100) ,HZ( 100) , VX( 100) , V¥( 100) ,VZ( 100) , DT( 100) ,NX(100) 
DIMENSION NY( 100) ,NZ( 100) ,RDV( 100) ,EX(100) ,EY( 100) ,EZ( 100) 
DIMENSION NDEC 100) ,EDR( 100) ,V( 100) ,HT( 100) ,RDRF( 100) , INTA( 100) 
DIMENSION RONA( 100) ,THJO( 100) , ATE( 100) ,CTE( 100) 
COMMON/OBLATE1/DAY ,RX,RY,RZ,VX,VY,VZ,DT,HX,HY,HZ,NX,NY,NZ,K,KK 
COMMON/OBLATE2/RDV,R,V,EX,EY,EZ2,MU,NDE,EDR,H,N,E,P,1,0N,W,F 
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aan 


20 


re 


COMMON/OBLATE3/PI,EC,M,A,HT,ER,TH,THD,RTD,MD,WD,OMD,ID,ECD 
COMMON/OBLATE4/FD, LINE, J,THF ,THFD,IF, IFD,OMF ,OMFD,RF,INT,ROM 
COMMON/OBLATE5/RJ,DR,DID,DTHD,DOMD , ESTERR, ACTERR , TERROR, JVER 
COMMON/OBLATE6/RFX,RFY,RFZ,ARC ,ARCD,RDRF ,DRV,RJ2,JN, JORBIT 
COMMON/OBLATE7 /INTA, ROMA, THJO, ATE , CTE 


CALCULATE R,V,E,H,N,P,1,0M,W,F,EC,M,HT,TH 


DO 10 J=1, KK 
CALL CROSS(RX(J) ,RY(J) ,RZ(J) , VX(J) , VYCS) ,VZ( J) ,HX(J) ,HY( J), 
HZ(J 


CALL CROSS(0. ODO,0. ODO, 1. ODO, HX(J) ,HY(J) ,HZ(J) ,NX(J) ,NY(J), 
NZ2(J)) 


CALL DOT(RX(J) ,RY(J) ,RZ(J) ,VX(J) ,V¥(J) ,VZ(T) ,RDV(J)) 
R(J)=DSORT(RX(J)*RX(I)+RY(J)*RY(J)4RZ(I)*RZ(J)) 
V(J)=DSORT( VX(J)*VX(J)FVY(J)#VY( J)4VZ(J)*VZ(J)) 
EX(J)=((V(J)*V(J)-MU/R(J) )*RX( J) -RDV(J)*VX(J))/NU 
EY(J)=((V(I)*V(J)-NU/R(J) )#RY(J) -RDV(J)*V¥()) /MU 
E2(J)=((V(3)*V(I)-MU/R( J) )*RZ(J) -RDV(T)*VZ(T)) /MU 
CALL DOT(NX(J) ,NY(J) ,NZ(J) ,EX(J) ,EY(J) ,EZ(J) ,NDE(J)) 
CALL DOT(EX(J) ,EY(J) ,EZ(J) ,RX(J) ,RY(J) ,RZ(J) ,EDR(J)) 
H(J)=DSQRT(HX(J)*HX(J)+HY(J)*HY(J)+HZ(3)*HZ(J)) 
N( J)=DSORT(NX(J)*NX(J)-4NY(J)*NY(J)-4NZ(J)*NZ(J)) 
E(J)=DSQRT(EX(J)*EX(J)+EY(J)*EY(J}+EZ(J)*EZ(J)) 
P(J)=H(J)*H( I) /MU 
1( J)=DACOS(HZ(J)/H(J)) 
OM(J)=DACOS(NX(J)/N(J)) 
W(J)=DACOS(NDE(J)/(N(J)*E(J))) 
F(J)=DACOS(EDR(J)/(E(J)*R(J)}) 
IF(NY(J). LE. 0. ODO) THEN 

OM(J)=2. 0DO*PI-OM(J) 
ENDIF 
IF(EZ(J). LE. 0. ODO) THEN 

W(J)=2. ODO*PI-W(J) 
ENDIF 
IF(RDV(J). LE. 0. ODO) THEN 

F(J)=2. ODO*PI-F(3) 
ENDIF 
EC( J)=DACOS( (E(J)+DCOS(F(J)))/(1. ODO+E(J)*DCOS(F(J)))) 
IF(F(J). GE. PI) THEN 

EC(J)=2. ODO*PI-EC(J) 
ENDIF 
M(J)=EC(J)-E(J)*DSIN(EC(J)) 
AC J)=P(J)/(1-E(J)*E(J)) 
AT=(MU/(N(1)*N(1)))***(1. 0D0/3. ODO) 
HT(J)=R(J)-ER | 
RJ=3. ODO#RI2“ER“ER/(2. ODO*P(1)*P(1)) 
TH(J)=F(J)+H(J) 
THD(J)=TH(J)*RTD 
IF(THD(J). GT. 360. ODO) THEN 

THD( J)=THD(J) -360. ODO 

GOTO 20 
END IF 
TH(J)=THD(J) /RTD 














G CONVERT ORBITAL ELEMENTS TO DEGREES 


G 
MDC J)=H(J)*RTD 
WDC J)=WCJ)*RTD 
OMD(J)=OM(J)*RTD 
ID(J)=ICJ)*RTD 
ECD( J)=EC(J)}*RID 
FD(J)=F(J)*RTD 
THD( J)=TH(J)*RID 
10 CONTINUE 
RETURN 
END 
C 
Cc kiviinvitedtciivciciicicicttea eR 
C ve * 
C * SUBROUTINE CROSS * 
Cc * * 
Cc dedevebieekionickiviicicdeiicteinice kets 
C 
SUBROUTINE CROSS(AX,AY,AZ,BX,BY,BZ,CX,CY,CZ) 
SMPLICIT DOUBLE PRECISION (A-I,M-Z) 
ea C 
C CALCULATE THE CROSS PRODUCT OF TWO VECTORS A AND B 
G 
CX=AY**BZ-AZ*BY 
CY=AZ*BX-AX*BZ 
CZ=AX*BY-AY*BX 
C 
RETURN 
END 
C 
Cc dntineivnlcleicteleiicicteiesseK se 
C ve " ve 
c te SUBROUTINE DOT t 
CG * vw 
C teieiedieteldekictsisticteleleicisiicieicictevckitetsss 
C 
SUBROUTINE DOT(AX, AY,AZ,BX,BY,BZ,ADB) 
IMPLICIT DOUBLE PRECISION (A-I,M-Z) 
G 
C CALCULATE THE DOT PRODUCT OF TWO VECTORS A AND B 
C 
ADB=AX**BX+AY*BY+AZ"BZ 
C 
RETURN 
END 
C 
Cc Selsieieleisieveiotetricisietcictvictviotcteiokicivisitstcts 
G Se we 
C : SUBROUTINE PINITIAL % 
Cc ve * 
Cc dedetededeisvetedetetetacetoteteisdsiciscetetssctesescterets seeds 
C 


SUBROUTINE PINITIAL 
IMPLICIT DOUBLE PRECISION (A-I,N-Z) 
CHARACTER*20 LINE 


110 











aQaag 


C 
C 
C 
C 
C 
C 
C 





beet 


tet t 





DIMENSION M(100) ,MD( 100) ,E( 100) ,W( 100) ,WD( 100) ,0OM( 100) ,OMD( 100) 
DIMENSION I(100),ID(100) ,F( 100) ,FD( 100) ,EC( 100) ,ECD( 100) , A( 100) 
DIMENSION R(100) ,H( 100) ,N(100) , TH( 100) , THD( 100) 

DIMENSION RF(100) ,IF(100) ,IFD( 100) ,OMF(100) ,OMFD( 100) , THF(100) 
DIMENSION THFD(100) ,P( 100) , JORBIT( 100) ,DR( 100) ,DID( 100) , DTHD( 100) 
DIMENSION DOMD( 100) ,RX( 100) ,RY( 100) ,RZ( 100) ,RFX(100) ,RFY( 100) 
DIMENSION RFZ(100),DRV(100) ,ARC( 100) , ARCD( 100) , DAY( 100) ,HX( 100) 
DIMENSION HY( 100) ,HZ( 100) ,VX( 100) , V¥(100) ,VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION NY(100) ,NZ( 100) ,RDV( 100) ,EX( 100) ,EY( 100) ,EZ( 100) 
DIMENSION NDE(100) ,EDR( 100) ,V(100) ,HT( 100) ,RDRF( 100) , INTA( 100) 
DIMENSION ROMA( 100) ,THJO( 100) ,ATE( 100) , CTE(100) 
COMMON/OBLATE1/DAY ,RX,RY,RZ, VX, VY,VZ,DT,HX,HY,HZ,NX,NY,NZ,K,KK 
COMMON/OBLATE2/RDV,R,V,EX,EY,EZ,MU,NDE,EDR,H,N,E,P,1,OM,W,F 
COMMON/OBLATE3/PI,EC,M,A,HT,ER, TH, THD,RTD,MD,WD,OMD,ID,ECD 
COMMON/OBLATE4/FD,LINE,J,THF,THFD,IF,IFD,OMF,OMFD,RF, INT, ROM 
COMMON/OBLATE5S/RJ,DR,DID,DTHD , DOMD ,ESTERR , ACTERR , TERROR , JVER 
COMMON/OBLATE6/RFX,RFY,RFZ,ARC,ARCD,RDRF ,DRV,RJ2,JN,JORBIT 
COMMON/ OBLATE7/INTA, ROMA, THJO, ATE , CTE 


PRINT INITIAL ORBITAL ELEMENTS 


WRITE(6,'(/)') 

WRITE(2,'(/)") 

WRITE(6,6000) ' ORBITAL ELEMENTS’ 
WRITE(2,6000) ' ORBITAL ELEMENTS' 
WRITE(6,6100) LINE 

WRITE(2,6100) LINE 


WRITE(6,6200) 'M ',MD(1),'N = ',N(1),'E = ',E(1), 





'y = ',WDC(1),'OM = ',OMD(1),'I = ',ID(1), 
‘EC = ',ECD(1),'A = ',A(1),'R = ',R(1), 
ans ( HDCA) =',HC),'F = ',FD(1), 
WRITE(2,6200) 'N = '.MD(1),'’N = ',N(1),"E = ',E(1), 
W = ‘WDC 1) "OM = ",OMD(1), "I = ',ID(1), 
1EC = ',ECD(1),'A = »AC1),°R = »R(1), 
pe = ar H = ,H(1), F = ,FD(1), 
= 3 
FORMAT( 3X, A) 
FORMAT(3X,A20/) 
FORMAT( 13(3X,A5,D18. 10/)/) 
RETURN 
END 
Seledeteleleledledeledesedededetedoiedededeldehehhk ented 
* te 
* SUBROUTINE INTEGRATE we 


* x I 


dededededededetesededesedece Tedd elle leLeLehe te Teleledehelehevede 


SUBROUTINE INTEGRATE 

IMPLICIT DOUBLE PRECISION (A-I,M-Z) 

CHARAGTER*20 LINE 

DIMENSION M( 100) ,MD( 100) ,E( 100) ,W( 100) ,WD( 100) ,0M( 100) ,OMD( 100) 
DIMENSION 1( 100) ,ID( 100) ,F( 100) ,FD( 100) ,EC( 100) ,ECD( 100) ,A( 100) 


lll 








DIMENSION R(100) ,H( 100) ,N( 100) ,TH( 100) , THD(100) 

DIMENSION RF(100) ,IF( 100) , IFD( 100) , OMF( 100) , OMFD( 100) , THF( 100) 
DIMENSION THFD(.*0) ,P( 100) , JORBIT( 100) ,DR( 100) ,DID( 100) ,DTHD( 100) 
DIMENSION DOME 0) ,RX( 100) ,RY( 100) ,RZ(100) ,RFX( 100) ,RFY(100) 
DIMENSION RFZ(? +5 ,DRV(100) ,ARC(100) , ARCD( 100) ,DAY( 100) ,HX(100) 
DIMENSION HY(10u) ,HZ( 100) , VX( 100) , V¥( 100) , VZ( 100) ,DT( 100) ,NX(100) 
DIMENSION NY( 100) ,NZ(100) ,RDV( 100) ,EX( 100) ,EY( 100) ,EZ(100) 
DIMENSION NDE(100) ,EDR( 100) ,V( 100) ,HT( 100) ,RDRF( 100) , INTA( 100) 
DIMENSION ROMA(100) , THJO( 100) , ATE( 100) , CTE( 100) 
COMMON/OBLATE1/DAY ,RX, RY RZ, VX, VY, VZ,DT, HX, HY ,HZ,NX,NY,NZ,K,KK 
COMMON/OBLATE2/RDV,R,V,EX,EY,EZ,MU,NDE,EDR,H,N,E,P,1,0M,W,F 
COMMON/OBLATE3/PI,EC,M,A,HT,ER, TH, THD,RTD,MD,WD,OMD, ID, ECD 
COMMON/OBLATE4/FD,LINE,J,THF ,THFD,IF,IFD,OMF,OMFD,RF, INT,ROM 
COMMON/OBLATES /RJ,DR,DID ,DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
COMMON/OBLATE6 /RFX,RFY ,RFZ,ARC,ARCD , RDRF ,DRV,RJ2, JN, JORBIT 
COMMON/OBLATE7 /INTA,ROMA,THJO, ATE, CTE 


EQUATE INITIAL VALUES 


aaa 


THE(1)=TH(1) 
THED( 1)=THD(1) 
IF(1)=1(1) 
IFD(1)=1D(1) 
OMF(1)=0M(1) 
OMFD( 1)=OMD(1) 
RF(1)=R(1) 


ESTIMATE UPPER BOUND OF THETA 
THF ( J)=TH(J)+(J-1)*0. 57D0*2. ODO*PI ROM00080 
INITIALIZE 


aQaaaqg aaa 


(1X10)-12 

ESTERR=0. 000000000001D0 
INT=DT(J)*H(1) 

DTHF=0. 1745329251994 

C (1X10)-10 

TERROR=0. 0000000001D0 


10 CALL ROMBERG ROM00190 
ACTERR=INT~-ROM 
IF( ACTERR. LT. 0. ODO) THEN 
THF(J)=THF( J) -DTHF 
GOTO 10 
ENDIF 
TEMPTHF=THF (J) 
GOTO 20 


30 CALL ROMBERG 
ACTERR=INT-ROM 
20 IF(ACTERR. GE. 0. ODO) THEN 
IF( (ACTERR/ INT). LE. ESTERR) THEN 
GOTO 40 
ELSE 
TEMPTHF=THF( J) 
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arAaggagaAaNa 


aaa a 


aaa 


‘ELSE 


THF ( J)=THF(J)+DTHF 
GOTO 30 
ENDIF 


DTHF=DTHF/2. ODO 
THF ( J)=TEMPTHF+DTHF 
GOTO 30 

END IF 


INTA( J)=INT 
ROMA( J)=ROM 


RETURN 

END 

PRT RRR EINER RIE ERE, 
* * 
* SUBROUTINE ROMBERG * 
* * 


WR IWATE TRIKE ERR CE ETE 


SUBROUTINE ROMBERG 

IMPLICIT DOUBLE PRECISION (A-I,N-2Z) 

CHARACTER*20 LINE 

DIMENSION M( 100) ,MD( 100) ,E( 100) ,WC 100) ,WD( 100) ,OMC 100) ,OMD( 100) 
DIMENSION I( 100) ,1D(100) ,F(100) ,FD(100) ,EC( 100) ,ECD( 100) ,A( 100) 
DIMENSION R( 100) ,H(100) ,N(100) ,TH¢ 100) , THD( 100) 

DIMENSION RF(100) ,IF( 100), IFD( 100) ,OMF(100) ,OMFD( 100) , THF( 100) 
DIMENSION THFD(100) ,P(100) , JORBIT(100) ,DR( 100) ,DID( 100) , DTHD( 100) 
DIMENSION DOND( 100) ,RX( 100) ,RY( 100) ,RZ(100) ,RFX(100) ,RFY( 100) 
DIMENSION RFZ(100) ,DRV( 100) ,ARC( 100) , ARCD(100) ,DAY( 100) ,HX( 100) 
DIMENSION HY( 100) ,H2( 100) , VX( 100) , VY( 100) ,VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION NY(109) ,NZ( 100) ,RDV( 100) ,EX(100) ,EY( 100) ,EZ( 100) 
DIMENSION NDE( 100) ,EDR(100).,V( 100) ,HT( 100) ,RDRF( 100) , INTA( 100) 
DIMENSION ROMA( 100) , THJO( 100) , ATE( 100) , CTE( 100) 
COMMON/OBLATE1/DAY,RX,RY,RZ,VX,VY,VZ,DT,HX,HY,HZ,NX,NY,NZ,K,KK 
COMMON/OBLATE2/RDV,R,V,EX,EY,EZ,MU,NDE,EDR,H,N,E,P,I,OM,W,F 
COMMON/OBLATE3/PI,EC,M,A,HT,ER,TH,THD,RTD,MD,WD,OMD,ID,ECD 
COMMON/OBLATE4/FD, LINE ,J,THF,THFD,IF, IFD,OMF ,OMFD,RF, INT,ROM 
COMMON/OBLATES /RJ,DR, DID, DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
COMMON/OBLATE6/RFX,RFY,RFZ,ARC,ARCD,RDRF , DRV,RJ2,JN, ORBIT 
COMMON/OBLATE? /INTA, ROMA, THJO, ATE ,CTE 


EXTERNAL FUNC 

INITIALIZE VARIABLES 

HS=THF( J) -THF(1) 
FUNCA=FUNC(RJ,A(1),I(€1),E(1) ,W(1) ,TH(1) ,TRF(1) ,JVER) 
FUNCB=FUNC(RJ,A(1),1(1),E(1),WC1),TH(1) , THE(J) , JVER) 
P(1)=HS*(FUNCA+FUNCB) /2. ODO 

JM=1 

BEGIN THE ROMBERG LOOP. 


DO 10 JN= 1, 100 
OLDINT=P(1) 








ROM00190 


ROM00460 
ROM00560 
ROM00570 
ROM00580 
ROM00590 
ROM00610 
ROM00620 


ROM00630 
ROM00640 








20 


30 


C 
C 
C 
C 
C 
C 
C 








HH=HS 
SS=0. ODO 
TT=THF(1)+HH/2. ODO 
DO 20 L=1, JM 
T=TT 
SS=SS+FUNC(RJ,AC1),1(1) ,E(1) ,W(1) ,TH(1) ,T, JVER) 
TI=TT+HH 
CONTINUE 
SUM=HH*SS 
P( JN+1)=(P( JN)+SUM) /2. ODO 
D=1. ODO 
DO 30 JK = JN, 1, -1 
D=4. ODO*D 
P( JK)=P(JK+1)+(P( JK+1) -P( JK) )/(D-1.-0D0) 
CONTINUE 
ERROR=( P(1)-OLDINT) 
IF(JN. GE. 10) THEN 
IF (DABS(ERROR/OLDINT). LE. TERROR) THEN 
GOTO 40 
ENDIF 
ENDIF 
HS=HS/2..0D0 
JM=IM2 
CONTINUE 
ROM=P( 1): 
RETURN 
END 


Inkiiokwbiicnekiciekiekkicek kK II 
¥ ca 
we FUNCTION FUNC vt 
v8 * 
detervederevededevedevedededededeetededededeven Kiclevedeseteteds 


FUNCTION FUNC(RJ,A1,11,E1,W1,TH1,THFJ,JVER) 
IMPLICIT DOUBLE PRECISION (A-1,0-2Z) 


EXTERNAL RADIUS 


S=DSIN(11) 
§2=8*§ 
$4=82*§2 
S6=54*S2 
E2=E1*E1 


RAD=RADIUS(RJ,A1,11,E1,W1,TH1,THFJ,JVER) 

F ( SOLUTION ) 

IF( JVER. EQ. 1) THEN 

Y1=112. 0D0-75. ODO*"S6+269. OD0*54-296, ODO*S2 
Y2=RI*THFI*( 2. 5DO*S2-2. ODO) 

Y3=2. ODO*W1-Y2 


Y4=24, 0DO*(5. ODO*S2-4. 0D0)*(5. ODO*S2-4. ODO) 
Y5=E2*52*( 14. ODO-15. OD0*S2)*( 15. ODO*S2-13. ODO) 


114 


ROM00670 
ROM00680 
ROM00690 
ROMO00720 
ROM00730 
ROM00740 
ROM00750 
ROM00760 
ROM00780 
ROMO00800 
ROM00830 
ROM00840° 
ROM00850 
RON00860 
ROM00870 
ROM00900 


ROM00930 


ROM00940 
ROM00950 
ROM00960 
ROM00970 


ROM00030 

















Y6=9. ODO*E2+34. ODO 
Y9=12. ODO*(5. ODO*S2-4. ODO) 
Y12=9. OD0*E2~-34. ODO 










t YF=(2. 5D0*S2-2. ODO)*(THFJ-TH1)+E2*Y1*DSIN( ¥2)*DCOS(Y3)/Y¥4 


C 
YS=Y5*DCOS( 2. ODO*W1)/(2. ODO*Y9)+ 
+ E1*$2*( 15. ODO*S2-13. ODO)*DCOS( TH1+W1)/2. ODO+ 
e + E1*S2*(15. OD0*S2-13. 0D0)*DCOS(3. ODO*TH1-W1)/6. ODO+ 
+ $2*(15. ODO*S2-13. ODO)*DCOS( 2. ODO*TH1) /2. ODO+ 
+ (5. ODO*Y6*S4+4. ODO*Y12*S2-56. ODO*E2) /96. ODO 





Y=THF J-W1+RI*YF+RI*RI*THFI*YS 





F=(2. OD0-3. OD0*S2)*DCOS( 2. ODO*THFI)/2. ODO+ 





















+ E1*(S$2-1)*DCOS(Y)+E1*(3. OD0-4. OD0*S2)*DCOS( Y+2. ODO*THF J) /6. ODO+ 
+  El*(1. 0D0-2. ODO*S2)*DCOS(Y-2. ODO*THFJ) /2. ODO+S2-1. ODO+ 
+  E2*§2%*(15. OD0*S2-14. ODO)*DSIN( Y2)*DSIN(Y3) /Y9+ 
+  $§2%*DCOS(2. ODO*TH1) /2. ODO+E1*S2*DCOS( 3. ODO*TH1-W1)/6. ODO+ 
+ £E1*S2*DCOS(TH1+W1)/2. 600 
C 
ENDIF 
C 
C F ( SIMPLIFIED SOLUTION ) 
C 
IFC JVER. EQ. 2) THEN 
Cc 
: F=(2. OD0-3. ODO0*S2)*DCOS( 2. ODO* THF) /2. ODO+ 
+ E1*(S2-1)*DCOS( THFJ-W1)+ 
+ E1*(3, ODO-4. ODO*S2)*DCOS(3. ODO*THFI-W1) /6. ODO+ 
+ E1%*(1. ODO-2. OD0*S2)*DCOS( THFI+W1) /2. ODO+S2-1. ODO+ 
? +  $§2*DCOS(2. ODO*TH1)/2. ODO+E1*S2*DCOS( 3. ODO*TH1-W1) /6. ODO+ 
+ E1*S2*DCOS(TH1+W1)/2. ODO 
C 
ENDIF 
C 
C F ( TWO BODY SOLUTION ) 
C 
IF( JVER. EQ. 3) THEN 
F=0. ODO 
ENDIF 
C 
C FUNCTION 
C 
FUNC=RAD*RAD*(1. ODO+RJ*F) 
C 
RETURN 
END 
C 
C dededededede dele lecele TISAI NLA NNER 
. Cc oi ve 
C ve FUNCTION RADIUS * 
CG ve Xe 
‘ Cc Yeledededesedeteleletelesesele Lele NI AATEC 
C 


FUNCTION RADIUS(RJ,A1,11,E1,W1,TH1,TFJ,JVER) 
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IMPLICIT DOUBLE PRECISION (A-1I,0-2Z) ROM00030 


aac 


CALCULATE E, SINE, AND COSINE FUNCTIONS 


S=DSIN(11) 

S2=S8*S 

S4=$2*§2 

56=84*§2 

C=DCOS(T1) 

C2=C*C 
SC=DSIN(1I1)*DCOS(I1) 
E2=E1*E1 

PO=A1*(1. OD0~E2) 


RADIUS BOTTOM ( SOLUTION ) 
IF(JVER. EQ. 1) THEN 


Q aan 


Y1=112. 0D0-75. 0D0*S6+260. OD0*S4-296. ODO*S2 
Y2=RI*TFI*(2.5D0*S2-2. ODO) 

Y3=2. ODO*W1-Y2 

Y4=24. ODO*(5. ODO*S2-4. 0D0)*(5. ODO*S2-4. GDO) 
Y5=E2*S2*( 14. ODO-15. OD0*S2)*( 15. ODO*S2-13. ODO) 
Y6=9. ODO*E2+34. ODO 

Y9=12. ODO*(5. ODO*S2-4. ODO) 


Y11=15. ODO*( 2. ODO+E2)*S4-14. ODO*( 4. ODO+E2)*S2+24. ODO 
Y12=9, ODO*E2-34. ODO 
YF=(2. 5D0*52-2. ODO)*(TFI-TH1)+E2*Y1*DSIN(¥2)*DCOS(Y3) /Y4 


a a a a 


YS=Y5*DCOS( 2. ODO*W1)/( 2. ODO*Y9)+ 
E1*S2*( 15. ODO*S2-13. OD0)*DCOS(TH1+W1) /2. ODO+ 
E1*52*(15. OD0*S2-13. 0D0)*DCOS(3. ODO*TH1-W1) /6. ODO+ 
$2*( 15. OD0*S2-13. ODO)*DCOS(2. ODO*TH1) /2. ODO+ 
(5. ODO*Y6*S44+4. ODO*Y12*52-56. ODO*E2) /96. ODO 


tet + 


Y=TFI-W1ItRI*YF+RI*RISTEINYS 


RF1=1. ODO-1. 5D0*S2+E2*(1. ODO-1. 25D0*S2) - 
((2. ODO+S. ODO*E2)*S52-2. ODO*E2)*DCOS( 2. ODO*TFJ) /12. ODO+ 
E2*(9. OD0*S2-8. OD0)*DCOS( 2. ODO*Y) /12. ODO+ 
E1*(-11. ODO*S2+6, 0D0)*DCOS( Y+2. ODO*TFJ)/24. ODO+ 
E2*(-3, ODO*S2+2. ODO)*DCOS( 2. ODO*Y+2. ODO*TFI) /24. ODO+ 
E2*( 3, ODO*S2-2. 0D0)*DCOS( 2. ODO*Y-2. ODO*TFJ) /8. ODO+ 
E1*Y11*DSIN(Y2)*DSIN( TFI+W1) /Y9 


thet t+ 


RF2=E2*S2*( 15. OD0*S2-14. 0D0)*DSIN( Y¥2)*DSIN(Y3)/( 0. 5DO*Y9) -~ 
E2*S2*DCOS(Y-TH1+3. ODO*W1)/16. ODO+ 
E2*(3. ODO*S2-2. 0D0)*DCOS(Y-3. ODO*TH1+3. ODO*W1)/24. DO- 
E2*S2*DCOS(Y-5. ODO*TH1+3. ODO*W1) /16. ODO+ 
E1*(3. ODO*S2-2. 0D0)*DCOS( Y-2. ODO*TH1+2. ODO*W1) /4. ODO- 
3. ODO*E1*S2*DCOS( Y-4. ODO*TH1+2. ODO*W1)/8. ODO- 
E1*(S2+1. 0D0)*DCOS(Y+2. ODO*W1)/4. ODO 


bee bet 
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RF3=((5. ODO*E2-2. 0D0)*S2-2. ODO*E2)*DCOS(Y+TH1+W1)/8. ODO+ 
((5. ODO*E2+6. 0D0)*S2-4. ODO*(E2+1. ODO) )*DCOS( Y-THI+W1) /4. ODO+ 
(2. ODO*E2-S2*(5. ODO*E2+14. ODO) )*DCOS(Y-3. ODO*THI+W1) /24. ODO+ 
E2*(9. ODO*S2-4. ODO)*DCOS( Y+3. ODO*TH1-W1) /48. DO+ 
E2*(6. ODO~-7. ODO*S2)*DCOS( Y+TH1-W1)/8. ODO+ 
E2*(4..0D0-5. 0DO*S2)*DCOS( Y-TH1-W1)/16. ODO+ 
E1*¢2. ODO*S2-1. ODO)*DCOS( Y+2. ODO*TH1) /4. ODO 


bette t+ 


RF4=E1*(1. OD0-3. ODO*S2)*DCOS(Y-2. ODO*TH1) /4. ODO+ 
E1*(2..0D0-3. 0D0*S2)*DCOS(Y) /4. ODO+ 
E1*S2*DCOS( TH1+W1)+S2*DCOS( 2. ODO*TH1)+ 
E1*S2*DCOS{ 3. ODO*TH1-W1)/3. ODO 


+e + 


RFB=1. ODO+E1*DCOS( Y)+RI*(RFI+RF2+RF3+RF4) 
ENDIF 

RADIUS BOTTOM ( SIMPLIFIED SOLUTION ) 

IF( JVER. EQ. 2) THEN 


RF1=1. ODO-1. 5D0*S2+E2*(1. ODO-1. 25D0*§2) - 
((2. ODO+5. ODO*E2)*S2-2. ODO*E2)*NCOS( 2. ODO*TFIJ) /12. ODO+ 
E2*(9. ODO*S2-8. O0D0)*DCOS( 2. ODO*(TFIJ-Wi))/12. ODO+ 
E1*(-11. OD0*S2+6. OD0)*DCOS(3. ODO*TFI-W1)/24. ODO+ 
E2*(-3. ODO*S2+2. ODO)*DCOS(4. ODO*TFI-2. ODO*W1)/24. ODO+ 
E2*(3. ODO*S2-2. OD0)*DCOS( 2. ODO*W1)/8. ODO- 
E2*S2*DCOS(TFJ-TH1+2. ODO*W1)/16. ODO 


ee 


RF2=E2*(3, ODO*S2-2. 0DO)*DCOS(TFI-3. ODO*TH1+2, ODO*W1)/24. DO- 
E2*S2*DCOS(TFI-5. ODO*TH1+2. ODO*W1)/16. ODO+ 
E1**(3. ODO*S2-2. ODO)*DCOS(TFI-2. ODO*TH14W1) /4. ODO- 
3, ODO*E1*S2*DCOS(TFJ-4. ODO#TH14W1) /8. ODO- 
E1**(S2+1. 0D0)*DCOS(TFI+W1) /4. ODO+ 
((5. ODO*E2-2. 0D0)*S2-2. ODO*E2)*DCOS( TFI+TH1) /8. ODO+ 
((5. ODO*E2+6. 0D0)*S2-4. ODO*(E2+1, ODO) )*DCOS(TFJ-TH1) /4. DO 


bette t 


RF3=( 2, ODO*E2-S2*(5. ODO*E2+14. ODO) )*DCOS( TFJ-3. DO*TH1) /24. DO+ 
E2*(9, OD0*S2-4. 0D0)*DCOS( TFJ+3. ODO*TH1-2. ODO*W1)/48. DO+ 
E2*(6. ODO~7. ODO*S2)*DCOS( TFJ+TH1-2. ODO*W1)/8. ODO+ 
E2*(4. OD0-5. ODO*S2)*DCOS( TFJ-TH1-2. ODO*W1)/16. ODO+ 
E1*( 2. ODO*52-1. ODO)*DCOS( TFJ+2. ODO*TH1-W1)/4. ODO+ 
Ei*(1. 0D0-3. 0D0*S2)*DCOS(TFJ-2. ODO*TH1-W1)/4. ODO+ 
E1*(2. ODO-3. ODO*S2)*DCOS( TFJ-W1) /4. ODO 


bet EH+ 


RF4=E1**S2*DCOS(TH1+W1)+52*DCOS( 2. ODO*TH1)+ 
+ E1*S2*DCOS( 3. ODO*TH1-W1)/3. ODO 


RFB=E 1*DCOS( TFJ-W1+RJ*( TFI-TH1)*( 2. 5D0*S2-2. ODO) )+ 
+ 1. ODO+RI*(RFAI+RF2+RF3+RFS) 


ENDIF 
RADIUS BOTTOM ( TWO BODY SOLUTION ) 


IF(JVER. EQ. 3) THEN 











RFB=1. ODO+E1*DCOS(TFJ-W1) 
ENDIF 


‘RADIUS 
RADIUS=P0/REB 


RETURN ROM00410 
END -ROMOO420 + 
RAR BT RR RAIA RRR 
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SUBROUTINE FORMULA 

IMPLICIT DOUBLE PRECISION (A-I,M-Z) 

CHARACTER*20 LINE ; 

DIMENSION M( 100) ;MD(100) ,E( 100) ,W( 100) ,WD( 160) ,OM( 100) , OMD( 100) 
DIMENSION 1I( 100) ,ID(100) ,F( 100) ,FD( 100) ,ECC 100) ,ECD( 100) ,A( 100) 
DIMENSION RC 100) ,H( 100) ,N( 100) , THC 100) , THD( 100) 

DIMENSION RF(100) ; I¥( 100) ,IFD( 100) ,OMF( 100) ,OMFD( 100) , THF( 100) 
DIMENSION THFD( 100) ,P( 100) , JORBIT( 100) ,DR( 100) ,DID( 100) ,DTHD( 100) 
DIMENSION DOMD( 100) ,RX( 100) ,RY( 100) ,RZ(100) ,RFX(100) ,RFY( 100) 
DIMENSION RFZ(100) ,DRV( 100) ,ARC( 100) , ARCD( 100) , DAY(.100) ,HX( 100) 
DIMENSION HY( 100) ,HZ( 100) , VX(100) , VY¥( 100) ,VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION NY¥( 100) ,NZ( 100)-.,RDV( 100) ,EX( 100) ,EY( 100) ,EZ2( 100) 2 
DIMENSION NDE( 100) ,EDR(100) ,V(100) ,HT( 100) ,RDRF(100) , INTA( 100) 
DIMENSION ROMA( 100) , THJO( 100) , ATE( 100) ,CTE( 100) 
COMMON/OBLATE1/DAY ,RX,RY,RZ,VX,VY,VZ,DT,HX,HY,HZ,NX,NY,NZ,K,KK 
COMMON/OBLATE2/RDV,R,V,EX,EY,EZ,MU,NDE,EDR,H,N,E,P,1,0M,W,F be 
COMMON/OBLATE3/PI ,EC,M,A,HT,ER,TH,THD,RTD,MD,WD,OMD,ID,ECD 
COMMON/OBLATE4/FD, LINE ,J,THF ,THFD,IF,IFD,OMF,OMFD,RF,INT,ROM 
COMMON/OBLATES /RJ,DR, DID, DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
COMMON/OBLATE6/RFX,RFY,RFZ,ARC,ARCD,RDRF ,DRV,RJ2,JN, JORBIT 
COMMON/OBLATE7/INTA, ROMA, THJO, ATE, CTE 


EXTERNAL RADIUS 

CALCULATE E, SINE, AND COSINE FUNCTIONS 

S=DSINCI(1)) 

S2=8*S 

$4=52*52 

S6=54"S2 

C=DCO8(I(1)) 

C2=C*C 

SC=DSIN(1I(1))*DCOS(I(1)) 

E2=E(1)*E(1) 

FORMULA ( SOLUTION ) =e 
IF(JVER. EQ. 1) THEN . 


Yi=112. OD0-75. ODO*S6+260. OD0*54~-296. ODO*S2 
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Y2=RI*THF(IJ)*( 2. 5D9*S2-2. ODO) 

Y3=2. ODO*W(1)-Y2 

Y4=24. 0DO0*(5. ODO*S2-4, ODO)*(5. ODO*S2-4. ODO) 

Y5=E2*S2**( 14. ODO-15. ODO*S2)*( 15, 0D0*S2-13. ODO) 

Y6=9. ODO*E2+34. ODO 

Y7=15. ODO*S4-45.-0D0%S2+28. ODO 

Y8=6. ODO*(5. ODO*S2-4. 0D0)*(5. ODO*S2-4. ODO) 4 
Y9=12. ODO*(5. ODO*S2~4. ODO) ; 


Y¥10=(6. ODO-S2)/12. ODO-E(1)*S2*DCOS(3. ODO*TH(1)-W(1))/3. ODO- 
+ S2*DCOS( 2. ODO*TH( 1) )+E2*(7. ODO*S2-4. ODO) /24. ODO 


Y¥12=9. ODO*E2-34. ODO 
YF=(2. 5D0*S2-2. ODO)*(THF( J) -TH( 1) )+E2*Y1*DSIN( Y2)*DCOS(Y3)/Y4 


YS=Y5*DCOS(2.'0D0*W(1))/(2..0D0*¥9)+ 
E(1)*S2*( 15. ODO*S2-13. ODO)*DCOS(TH(1)+W(1))/2. ODO+ 
E(1)*S2*(15. ODO*S2-13. ODO)*DCOS(3. ODO*TH( 1)-W(1))/6. ODO+ 
$2**( 15. ODO*S2-13. ODO)*DCOS( 2. ODO*TH(1))/2. ODO+ 
(5. ODO*Y6*S444. ODO*Y12*§2-56. ODO*EZ) /96. ODO 


bb + 


Y=THF (J) -WC1)+RISYF+RI*RI*THF( J)*YS 
CALCULATE INCLINATION ( SOLUTION ) 
IF(J)=1(1)+SC*RI*(DCOS( 2. ODO*THF(J))/2. ODO+ 

E2*(14. DO-15. D0*S2)*DSIN( Y2)*DSIN(Y3) /(12. DO*(5. DO*S2-4. DO) )- 


DCOS( 2. ODO*TH(1))/2. ODO-E( 1)*DCOS( 3. ODO*TH( 1)-W(1))/6. ODO- 
E(1)*DCOS( THC 1)+W(1))/2. ODO) 


bet + 


CALCULATE LONGITUDE OF THE ASCENDING NODE ( SOLUTION ) 
JQ=RISRI 


OMF( J)=OM( 1)+C*RJ*( TH( 1) -THF(J)+DSIN( 2. ODO*THF(J))/2. ODO- 
E(1)*DSINCY)+E(1)*DSIN(Y+2. ODO*THF( J) )/6. ODO- 
E(1)*DSINCY-2. ODO*THF(J))/2. ODO-DSIN( 2. ODO*TH(1))/2. ODO+ 
E(1)*DSIN(TH(1)-W(1))-E( 1)*DSINC3. ODO*TH( 1)-W(1))/6. ODO- 
E(1)*DSIN(TH(1)+W(1))/2. ODO+E2*¥7*DSIN(Y2)*DCOS(Y3) /Y8)+ 
C*RI2*THF( J)*(E2*S2*( 15. ODO*S2-14. ODO)*DCOS( 2. DO*W(1))/Y9- 
E(1)*S2*DCOS(TH(1)+W(1))+Y10) 


bE EH 


ENDIF 

FORMULA ( SIMPLIFIED SOLUTION ) 

IF( IVER. EQ. 2) THEN 

CALCULATE INCLINATION ( SIMPLIFIED SOLUTION ) 
IF(J)=1(1)+SC*RI*(DCOS( 2. ODO*THF( J) )/2. ODO+ 


+ E(1)*DCOS(3. ODO*THF(J)-W(1))/6. ODO+ 
+ E(1)*DCOS( THF(J)+W(1))/2. ODO-DCOS( 2. ODO*TH(1))/2. ODO- 
+ E(1)*DCOS(3. ODO*TH(1)-W(1))/6. ODO- 
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E(¢1)*DCOS(TH(1)+W(1))/2. ODO) 


CALCULATE LONGITUDE OF ‘THE ASCENDING NODE ( SIMPLIFIED SOLUTION ) 

OMF(3)=OM(1)+C*PJ**( THCL)-THF(J)+DSIN( 2. ODO*THF(J))/2. ODO- 
EC 1)*DSINCTHE( J) -W(1))+E(1)*DSIN(3. ODO*THF(J)-W(1))/6. ODO+ 
E(1)*DSIN(THF(J)+W(1))/2Z. ODO-DSIN(2. ODO*TH(1))/2. ODO+ 
E(1)*DSIN(TH( 1)-¥(1))-E(1)*DSIN(3. ODO*TH(1)-W(1))/6. ODO- 
E(1)*DSIN(TH(1)+W(1))/2. ODO) 

ENDIF 

FORMULA ( TWO BODY SOLUTION ) 

IF(JVER. EQ. 3) THEN 

CALCULATE INCLINATION ( TWO BODY SOLUTION ) 

IF(J)=1(1) 


CALCULATE LONGITUDE OF THE ASCENDING NODE ( TWO BODY SOLUTION ) 


-OMF(J)=0M(1) 


ENDIF 

CALCULATE RADIUS ( SOLUTION, SIMPLIFIED, OR TWO BODY ) 
RFCJ)=RADIUS(RJ,AC1), 101) ,E(1) ,WC1), THC1) , THF( J) , JVER) 
CONVERT ANGLES TO DEGREES 


OMED( J)=OMF(J)*RTD 

IFD(J)=1F(J)*RTD 

THED(J)=THF(J)*RTL 

JORBIT(J)=0 

IF(THFD(J). GT. 360. ODO) THEN 
THED(J)=TEFD(.7) 360. ODO 
JORBIT(J)=JORBIT(I)+1 
GOTO 10 

END IF 

THJO(J)=JORBIT(J)*2. ODO*PI+TH( J) -TH(1) 

IF(OMFD(J). GT. 360. ODO) THEN 
OMFD( J)=OMFD(J)~360, ODO 
GOTO 20 

END IF 

THF(J)=THFD( J) /RTD 

OME(J)=OMFD( J) /RTD 


CALCULATE DELTAS 


DR(J)=RF(J)-RCJ) 

DID(J)=IFD(J)-ID(J) 

DTHD( J)=THFD( J) -THD(J) 

IF(DABS(DTHD( J) ). GE. 180. OD0)THEN 
IF(DTHD( J). LT, 0. ODO)THEN 
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DTHD( J)=DTHD( J)+360. ODO 
ELSE 


_ DTHD(J)=DTHD(J)-360. ODO 
ENDIF 

ENDIF 

DOMD( J)=OMFD(J)-OMD( J) 


RETURN 
END 


KIMANI REAR INARA ERE 
* * 
SUBROUTINE INERTIAL * 
te + 
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SUBROUTINE INERTIAL 

IMPLICIT DOUBLE FRECISION (A-I,M-Z) 

CHARACTER*20 LINE : 
DIMENSION M(100) ,MDC100) ,E( 100) ,W(100) ,WD( 100) ,OM( 100) , OMD( 100) 
DIMENSION I(100) ,ID( 100) ,F(100), FD(100), EC(100), ECD(100), A(100) 
DIMENSION R( 100) ,H(100).; -N( 100), TH( 100) , THD( 100) 

DIMENSION RF( 100) , IF(100) , IFD( 100), OMF( 100) , OMFD( 100) , THF( 100) 
DIMENSION THFD( 100) , P(100), JORBIT( 100), DR( 100), DID( 106) , DTHD( 100) 
DIMENSION DOMD( 100) ;RX( 100) ,RY(100) ,RZ( 100) ,REX( 100) ,RF¥( 100) 
DIMENSION RFZ( 100) ,DRV(100) , ARC( 100) , ARCD( 100) ,DAY( 100) , HX( 100) 
DIMENSION HY(100), HZ( 100), ¥X(100) , v¥(100), V2( 100), DT( 100), NX( 100) 
DINENSION NY(100) ,NZ(100), RDV(100), EX(100), EY(100), EZ(100) 
DIMENSION NDE(100), EDR(100) , V(100) ,HT( 100), RDRF(100), INTA( 100) 
DIMENSION EARC( 100) , EARCD( 100), PDR( 100) 

DIMENSION ROMA( 100) , THJO( 100) , ATE( 100) ,CTE( 100) 
COMMON/OBLATE1/DAY ,RX¥,RY,RZ,VX,VY,VZ,DT,HX,HY,HZ,NX,NY,NZ,K,KK 
COMMON/OBLATE2/RDV,R,V,EX,EY,EZ2,MU,NDE,EDR,H,N,E,P,1,OM,W,F 
COMMON/OBLATE3/PI,EC,M,A4,HT,ER,TH,THD,RTD,MD,WD,OMD,ID,ECD 
COMMON/OBLATE4/FD,LINE,J,THF,THFD, IF ,IFD,OMF,OMFD,RF,INT, ROM 
COMMON/OBLATE5/RJ,DR,DID,DTHD ,DOMD, ESTERR , ACTERR, TERROR , JVER 
COMMON/OBLATE6/RFX,RFY,RFZ,ARC,ARCD,RDRF,DRV,RJ2,JN,JORBIT 
COMMON/OBLATE7 /INTA, ROMA, THJO, ATE ,CTE 
COMMON/SPECIAL/EARC , EARCD , PDR, ENG, ENGF 


CALCULATE INITIAL ENERGY 


ENG=V(1)*V(1)/2. ODO-MU/R(1) 
ENGF=V(1)*V(1)/2. ODO-MU/RF(1) 


CALCULATE INERTIAL COORDINATES 
RFX(J)=RF(J)*(DCOS( THF( J) )*DCOS(OMF(J))- 
DSINC(THF( J) )*DCOS( IF( J) )*DSIN(OMF(J))) 
RFY(J)=RF(J)**(DCOS( THF( J) )*DSIN( OMF( J) )+ 
DSIN( THF( J) )*DCOS( IF( J) )*DCOS( OMF( J) )) 
RFZ(J)=RF(J)*(DSIN( THF(J))*DSINCIF(J))) 
CALCULATE DR VECTOR 


DRV(J)=DSQRT( (RFX(J) -RX(J) )*(RFX(J) -RX(J) + 
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+ (RFY(J).-RY(J) )*(RFY(J)-RY(J))+ 
+ (RFZ(J)-RZ(J) )*(RFZ(J)-R2(J))) 
PDR(J)=DRV(J)/R(J) 


CALCULATE ANGLE ERROR 


aaqga 


‘CALL DOT(RX(J) ,RY(J) ,RZ(J) ,REX(J) ,RFY(J) ,RFZ(J) ,RDRF(J)) 
ARC( J)=DACOS(RDRF(J)/(R(J)*RF(J))) 

CC=RF(J) 

CCP=R(J) 

BB=ER 

AA=DSQRT(BB*BB+CC*CC-2. 0DO*BB*CC*DCOS(ARC(J)/2. GDO)) 
AAP=DSQRT( BB*BB+CCP*CCP-2. ODO*BB*CCP*DCOS(ARC(J) /2. 0D0)). 
CCA=PI-DASIN(CC*DSIN(ARC(J)/2. 0D0)/AA) 
CCPA=PI-DASIN(CCP*DSIN(ARC( J) /2. ODO) /AAP) 

EARC(J)=2. ODO*PI-CCA-CCPA 

ARCD(J)=ARC(J)*RTD 

EARCD(J)=EARC(J)*RTD 


CALCULATE DOWNRANGE AND CROSSRANGE ERRORS 


aaa 


ATE( J)=R(J)*(DTHD(J) /RID+DCOS( I( J) )*DOMD(J) /RTD) 
CTE( J)=R(J)*(DSINCTH( J) )*DID( J) /RTID- 
+ DCOS(TH( J) )*DSIN( IC J) )*DOMD( J) /RTD) 


RETURN 
END 
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SUBROUTINE RESULTS 

IMPLICIT DOUBLE PRECISION (A-I,M-Z) 

CHARACTER*20 LINE 

CHARACTER*11 VERSION 

DIMENSION M(100) ,MDC100) ,E( 100) ,W(100) ,WD( 100) ,OM( 100) , OMD( 100) 
DIMENSION I( 100) ,ID( 100) ,F( 100) ,FD( 100) ,EC( 100) ,ECD( 100) ,A( 100) 
DIMENSION R( 100) ,H( 100) ,N( 100) , TH( 100) , THD( 100) 

DIMENSION RF(100) , IF( 100) ,IFD( 100) ,OMF( 100) ,OMFD( 100) , THF( 106) 
DIMENSION THFD( 100) ,P( 100) , JORBIT( 100) ,DR( 100) ,DID( 100) , DTHD( 100) 
DIMENSION DOMD( 100) ,RX( 100) ,RY( 100) ,RZ( 100) ,RFX( 100) ,RFY( 100) 
DIMENSION RFZ(100) ,DRV(100) ,ARC( 100) , ARCD( 100) ,DAY( 100) ,HX( 100) 
DIMENSION HY( 100) ,HZ( 100) , VX( 100) , V¥( 100) ,VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION NY(100) ,NZ( 100) ,RDV( 100) ,EX( 100) ,EY( 100) ,E2( 100) 
DINENSION NDE( 100) ,EDR( 100) ,V( 100) ,HT( 100) ,RDRF( 100) , INTA( 100) 
DIMENSION EARC( 100) ,EARCD(100) , PDR( 100). 

DIMENSION ROMA( 100) , THJO( 100) , ATE( 100) ,CTE( 100) 
COMMON/OBLATE1/DAY ,RX,RY,RZ,VX,VY,VZ,DT,HX,HY ,HZ,NX,NY,NZ,K,KK 
COMMON/OBLATE2/RDV,R,V,EX,EY,EZ,MU,NDE,EDR,H,N,E,P,1,OM,W,F 
COMNON/OBLATE3/PI,EC,M,A,HT,ER,TH,THD,RTD,MD,WD,OMD,ID,ECD 
COMMON/OBLATE4/FD, LINE ,J,THF ,THFD,IF,IFD,OMF,OMFD,RF,INT,ROM 
COMMON/OBLATES /RJ,DR,DID, DTHD, DOMD, ESTERR, ACTERR, TERROR, JVER 
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COMMON/OBLATE6/RFX,RFY,RFZ,ARC ,ARCD,RDRF ,DRV,RJ2, JIN, JORBIT 
COMMON/OBLATE7/INTA, ROMA, THJO, ATE , CTE 
COMMON/SPECIAL/EARC , EARCD , PDR, ENG , ENGF 


DR(1)=0. ODO 
DID(1)=0. ODO 
DTHD(1)=0. ODO 
DOMD(1)=0. ODO 
» DRV(1)=0. ODO 
ARCD(1)=0. ODO 
EARCD( 1)=0. ODO 
PDR(1)=0. ODO 
THJO( 1)=0. ODO 
ATE(1)=0. ODO 
CTE(1)=0. ODO 





IF(JVER: EQ. 1) THEN 
VERSION=' SOLUTION : 

ELSEIF( JVER. EQ. 2) THEN 
VERSION=' SIMPLIFIED’ 

ELSEIF( JVER. EQ. 3) THEN 
VERSION=' SECULAR ' 

ENDIF 

IF(RJ. EQ. 0. ODO)THEN 
VERSION='TWO BODY ' 

ENDIF 


OUTPUT RESULTS FOR DISSPLA 


aaa 





IF(JVER. EQ. 1) THEN | 

WRITE(3,3000) K | 

1 WRITE(3,3100) RJ | 
ENDIF 


pO 10 J = 1, KK 
WRITE(3,3100) DR(J),DID(J) ,DTHD(J) 
WRITE(3,3100) DOMD(J) ,DRV(J) ,EARCD(J) 
WRITE(3,3100) PDR(J) ,ATE(J) ,CTE(J) 
: WRITE(4,3100) THJO(J) 
0 CONTINUE 


PRINT RESULTS 


aQaare 


WRITE(6,'(/)') 
WRITE(2,'(/)') 
WRITE(6,6000) ' RESULTS‘ 
WRITE(2,6000) ' RESULTS' 
WRITE(6,6100) LINE 
WRITE(2,6100) LINE 
WRITE(6,6200) 'J = ',RJ 
7 WRITE(2,6200) ‘J = ',RJ 


WRITE(6,6300) ‘VERSION = ', VERSION 


WRITE(2,6300) ‘VERSION = ', VERSION 
P WRITE(6,6100) LINE 
WRITE(2,6100) LINE 











DO 20 J = K, KK 


WRITE(6,6400) ‘POINT = ',J,'ORBIT = ',JORBIT(J), 

*ROMBERG ITERATIONS = ',JN 
WRITE(2,6400) ‘POINT = ',J,'ORBIT = ',JORBIT(J), 

"ROMBERG ITERATIONS = ',JN 
WRITE(6,6500) 'R = ',R(J),'RF = ',RF(J),'DR = ',DR(J) 
WRITE(2,6500) 'R = ',R(J),'RF = ',RF(J),'DR = ',DR(J) 
WRITE(6,6500) 'I = ',ID(J),‘IF = ',IFD(J),'DI = ',DID(J) 
WRITE(2,6500) 'I = ',ID(J),‘IF = ',IFD(J),'DI = ',DID(J) 
WRITE(6,6500) ‘TH = ' ;THD( J), ‘THF = ' ;THED(J) ,'DTH = ' ,DTHD(J) 
WRITE(2,6500) "TH = °,THD(J),'THF = ',THFD(J),’DTH = ',DTHD(J) 
WRITE(6,6500) "OM = ',OMD(J),'OMF = ',OMFD(J),'DOM = ' ,DOMD(J* 
WRITE(2,6500) 'OM = ',OMD(J),‘OMF = ',OMFD(J),'DOM = * ,DOMD(J} 

t —t t Poe | t ae | 
WRIe,e800) EMCI = MCRL = TRC 

’ rT 3 )s -- sRYCJ), RZ _ »RZC(J) 

WRITE(6,6500) "RFX = ',RFX(J),'RFY = ',RFY(J),'RFZ = ',RF2(J) 
WRITE(2,6500) ‘RFX = *,REX(J),'RFY = ',RFY(J),'RFZ = ',RFZ(J) 
WRITE(6,6500) 'DRV = ',DRV(J),'PDR = *,PDR(J), 

"EARC= ' ,EARCD(J) 
WRITE(2,6500) ‘DRV = ',DRV(J),'PDR = ',PDR(J), 

*EARC= ‘' ,EARCD(J) 
WRITE(6,6500) ‘RTE = *,DR(J),'ATE = ',ATE(J), 

"CTE = ' ,CTE(J) 
WRITE(2,6500) ‘RTE =  DR(J), ‘ATE = ',ATE(J), 

"CTE = * ,CTE(J) 
WRITE(6,6600) ‘INT = ',INTA(J),‘ROM = ",ROMA(J) 
WRITE(2,6600) "INT = ',INTA(J),'ROM = ',ROMA(J) 

CONTINUE 

WRITE(6,6500) 'EG = ',ENG,'EGF = ',ENGF 
WRITE(2,6500) ‘EG = ',ENG,'EGF = ',ENGF 


WRITE(6,'(/)') 
WRITE(2,'(/)') 


FORMAT( 3X, £3) 
FORMAT(3(3X,D18. 10)) 


FORMAT(3X,A) 

FORMAT(3X,A20//) 
FORMAT(3X,A,F8. 6) 
FORMAT(3X,A,A11//) 
FORMAT(2( 3X, 48, 13/) ,3X,A21,13//) 
FORMAT(3(3X,A6 ,F23. 15/)) 
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6600 FORMAT(3(3X,A6,F23. 8/)) 
C 


RETURN 
END 
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